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Materials  Design  Principles  for  the 
Dynamic  Fracture  of  Laminar  Composite  Structures 

Final  Progress  Report  for  the  Period  August  15,  2005  -  February  28,  2006 

1.  Foreword 

Crack  bridging  (e.g.,  from  stitches  or  pins)  and  friction  have  profound  and  potentially 
very  useful  effects  on  delamination  crack  growth,  controlling  growth  rates  (damage 
levels)  and  the  energy  absorbed.  However,  the  implications  for  structural  design 
principles  remain  quite  obscure.  The  difficulty  is  that  no  simple  analogue  of  crack 
toughness,  which  underpins  static  structural  design,  exists  for  dynamic  cases  with  large 
scale  bridging  effects.  The  external  shape  of  the  structure  and  the  loading  configuration 
dictate  stress  waves,  frictional  contact  zones,  and  crack  tip  stress  intensity  factors  in  a 
way  that  is  very  difficult  to  approach,  other  than  by  brute,  case  specific  numerical 
simulation.  The  problem  is  compounded  by  the  reality  of  multiple  cracking,  a 
complexity  that  is  rarely  entertained  in  laboratory  fracture  specimen  design.  Physically 
sound  materials  models  for  the  important  structural  problem  of  multiple,  nonlinear 
cracking  in  laminated  structures  with  large-scale  bridging  due  to  friction  and 
reinforcement  remain,  in  spite  of  the  technological  importance  of  these  systems, 
undeveloped. 

We  are  conducting  a  program  of  basic  research  to  develop  engineering  principles  for 
dealing  with  dynamic,  multiple  cracking  damage  in  laminated  structures,  including  large 
scale  crack  bridging,  due  to  through-thickness  reinforcement,  and  friction.  Bridging  and 
friction  will  be  treated  by  materials  models  at  the  smallest  scales  relevant  to  the 
mechanisms.  By  reference  to  the  fundamentals  of  the  dynamic  growth  of  single  cracks, 
which  is  already  largely  understood  (although  some  interesting  mysteries  remain),  simple 
approaches  will  be  formulated  to  calculating  the  development  of  distributed  delamination 
cracks  in  laminated  structures  with  non-trivial  geometry  and  general  loading  conditions. 
To  treat  large  scale  bridging  effects,  structural  sub-component  models  must  support 
dimensions  of  ~  100  mm  or  more.  Our  approach  must  bridge  scales  ranging  from  this 
characteristic  size  down  to  that  of  micromechanisms  (friction,  fiber  bridging)  within  the 
process  zone  of  a  single  crack.  By  doing  this,  a  direct  link  will  be  established  between 
structural  performance  and  materials  design. 

Our  research  aims  to  create  a  systematic  method  for  simplified  design  of  laminated 
engineering  structures,  which  will  impact  the  design  and  performance  of  all  lightweight 
military  vehicles  and  structures.  Newly  gained  understanding  will  point  to  significant 
improvements  in  impact  and  ballistic  resistance  via  materials  and  structural  design. 

The  program  is  being  conducted  with  formal  collaborations  with  universities  in  the  U.S. 
and  Italy,  with  arrangements  pending  with  universities  in  Australia  and  Denmark.  All 
these  collaborations  involve  Ph.D.  candidate  students,  including  extended  stays  for 
students  at  Rockwell  Scientific. 
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Figure  1.  (a)  X-ray  image  of  delamination  crack  growing  from  the  tip  of  a 
sharp  notch  (Spearing  and  Beaumont,  1992).  (b)  Prediction  in  which 

delamination  and  splitting  cracks  are  the  only  damage,  (c)  Prediction  in 
which  delamination  and  splitting  cracks  are  accompanied  by  distributed 
softening  (due  to  microcracking)  in  the  laminate  plies. 

Figure  2.  Normalized  energy  release  rate  of  a  single  centrally  located 
stationary  crack  in  a  clamped  beam  specimen,  with  first  natural  period  T\ 
subject  to  different  pulse  durations,  Tp. 
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4.  Statement  of  the  Problem  Studied 

In  this  three-year  program  of  basic  research  we  shall: 

(i)  formulate  and  solve  multiple  delamination  problems,  by  methods  applicable  to 
reasonably  complex  structures  and  loading,  that  correctly  describe  fundamental 
aspects  of  dynamic  delamination  crack  growth  in  the  presence  of  large  scale 
bridging  and  friction; 

(ii)  establish  engineering  principles  for  modeling  and  designing  laminated  structures 
with  energy  absorption  optimized  by  tailoring  crack  bridging  and  friction 
mechanisms. 

The  main  thrust  of  our  work  will  be  theoretical.  Use  will  be  made  where  possible  of  data 
in  the  literature. 

Our  objectives  are  to: 

Conduct  basic  research  into  the  development  of  distributed  dynamic  delamination 
cracks  in  laminated  structures  with  non-trivial  geometry  and  general  loading 
conditions  and  in  the  presence  of  friction  and  bridging  due  to  through-thickness 
reinforcement  such  as  stitches  or  rods. 

Develop  engineering  principles  for  dynamic,  multiple  cracking  damage  in  laminated 
structures,  including  large  scale  crack  bridging  and  friction. 

Create  a  systematic  method  for  simplified  design  of  laminated  engineering  structures 
containing  through-thickness  reinforcement,  which  will  shorten  the  design  cycle  for 
all  lightweight  military  vehicles  and  structures. 

Indicate  means  for  achieving  significant  improvements  in  impact  and  ballistic  resistance 
via  materials  and  structural  design,  especially  via  manipulation  of  friction  or  the  bridging 
effects  of  through-thickness  reinforcement. 
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5.  Summary  of  the  Most  Important  Results 

Our  objective  in  the  first  year  is  to  formulate  multiple  delamination  problems  in  a  number 
of  cases  that  will  be  widely  representative  of  common  structural  applications. 
Formulations  are  sought,  using  approximations  where  they  are  valid  and  helpful,  that 
decouple  extrinsic,  shape  and  load  dependent  factors  from  the  intrinsic  nonlinear  material 
response  of  the  damaged  laminate.  Large  scale  bridging  conditions  of  friction  and 
through-thickness  reinforcement  will  be  included. 

Three  modeling  approaches  are  of  possible  interest,  viz.,  1)  beam  or  plate  level 
approximations;  2)  weight  function  methods;  and  3)  computational  methods.  In  the  first 
year,  selected  problems  will  being  solved  in  the  static  limit,  building  on  our  prior  related 
work  (Andrews  et  al.,  2006;  Yang  and  Cox,  2005).  Dynamic  solutions  are  also  being 
formulated  and  solved. 

Mixed  Mode  Cohesive  Models  -  Interacting  Failure  Mechanisms 

In  work  performed  under  other  programs,  Qingda  Yang  and  Brian  Cox  have  formulated 
high  fidelity  simulations  of  multiple  delamination  events  under  static  loading  in 
laminated  composites,  using  cohesive  zone  representations  of  the  nonlinear  processes 
associated  with  material  failure  in  the  crack  tip  region.  While  a  vast  and  long  literature 
exists  on  delamination  mechanics  (see  our  full  review  in  the  Introduction  to  (Yang  and 
Cox,  2005)),  our  work  presents  important  new  advances  that  make  a  big  difference  to 
whether  the  outcome  will  be  truly  high  in  fidelity.  In  our  new  ARO  program,  we  have 
extended  our  formulation  to  include  multiple  mechanisms  of  damage  and  shown  that  the 
interaction  between  different  mechanisms  can  dominate  damage  evolution. 

Key  points  we  have  contributed,  including  prior  work,  which  we  list  here  for 
completeness  in  presenting  our  point  of  view,  and  work  done  in  our  new  program,  are  as 
follows. 

1.  Cohesive  models  of  dominant  cracks  give  a  physically  consistent  description  of 
damage  initiation  (requiring  no  assumed  initial  defect),  the  progression  of  damage  to  a 
traction-free  crack,  and  the  propagation  of  the  crack  towards  the  limit  of  linear  elastic 
fracture  mechanics.  All  of  these  stages  of  crack  development  are  contained  in  a  single, 
unifying  model  -  the  cohesive  traction  law. 

2.  Linear  elastic  fracture  mechanics  (LEFM),  e.g.,  as  embodied  in  the  virtual  crack 
closure  technique,  is  very  limited  in  its  ability  to  generate  high  fidelity  simulations  of 
damage  evolution.  There  remains,  in  our  opinion,  no  rational  method  of  predicting 
damage  initiation  at  locations  of  elastic  stress  singularity  using  LEFM;  one  always  runs 
into  the  vexed  questions  of  what  initial  defect  might  exist  and  how  stress  singularities  are 
to  be  treated  computationally.  These  difficulties  are,  in  principle,  completely  overcome 
in  a  cohesive  zone  model,  which  provides  a  physically  valid  depiction  of  nonlinear 
material  behavior  prior  to  crack  formation;  and  introduces  a  length  scale,  the  cohesive 
zone  length,  which  plays  an  essential  role  in  developing  physically  appropriate, 
nonsingular  stress  representations. 
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3.  The  initiation  of  damage  and  the  evolution  of  the  crack  are  complex  processes, 
involving  continuous  changes  in  the  shape  of  the  crack  and  the  mode  ratio,  which  varies 
strongly  around  the  crack  front.  A  high-fidelity  simulation  must  be  allowed  to  predict  the 
crack  shape  evolution  -  it  cannot  be  specified  a  priori.  Shape  and  mode  ratio 
characteristics  follow  as  computational  results,  computed  in  real  time  during  a 
simulation,  naturally  and  very  easily  from  a  mixed  mode  cohesive  zone  representation. 

4.  The  cohesive  zone  has  a  characteristic  length  that  can  be  predicted  a  priori  using 
analytical  results,  which  we  and  our  collaborators  have  previously  derived  (Cox  and 
Marshall,  1994;  Massabo  and  Cox,  1999).  We  have  now  assembled  the  prior  results  into 
simple  rules  for  establishing  the  maximum  computational  mesh  size  for  simulations 
(Yang  and  Cox,  2005).  If  the  mesh  is  too  coarse  (elements  larger  than  the  cohesive  zone 
length),  mesh-independence  cannot  be  achieved.  Almost  all  prior  literature  on  cohesive 
zones  in  delamination  based  on  plate  elements  (or  similar  2D  elements)  fails  to  meet  this 
criterion!  Fortunately,  the  progression  of  computational  power  makes  it  feasible  to  use 
elements  that  are  properly  sized;  the  motivation  for  using  plate  elements  is  not  as 
compelling  as  it  has  been  in  the  past. 

5.  For  components  of  complex  shape  (e.g.,  containing  cut-outs),  strong  through¬ 
thickness  stresses  provide  a  further  reason  for  avoiding  plate  elements.  Three- 
dimensional  solid  finite  elements  are  preferable.  Once  again,  advances  in  computer 
power  make  using  solid  elements  feasible. 

5.  (New  work  in  this  program)  The  interaction  between  different  damage  modes  can 
be  very  strong.  High  fidelity  simulations  require  the  treatment  of  dominant  delamination 
cracks,  distributed  matrix  microcracking,  and  matrix  shear  damage  all  at  once. 

The  last  point  summarized  above  was  demonstrated  in  some  of  the  first  work  done  in  this 
ARO  program.  Details  are  shown  in  Figure  1.  In  earlier  work  of  ours  in  which  only 
delamination  and  splitting  cracks  were  included  as  possible  failure  modes,  many  of  the 
features  of  a  delamination  crack  emanating  from  a  sharp  slit  were  qualitatively 
reproduced,  but  the  exact  predicted  shape  of  the  delamination  crack  (Figure  lb)  was  not 
the  same  as  seen  experimentally  (Figure  la).  When  a  continuum  damage  representation 
of  ply  softening  was  included  along  with  interlaminar  and  splitting  cracks,  the  agreement 
with  experiment  became  excellent  (compare  Figures  la  and  lc).  These  results,  we 
believe,  strongly  encourage  the  prospect  of  high  fidelity  simulations  of  damage  using 
recent  physical  formulations  of  damage  and  current  computer  power. 

One  of  our  objectives  in  the  new  program  is  to  develop  similar  simulations  to  those 
depicted  in  Fig.  1  in  the  dynamic  regime. 

Beam  Theory  Solutions  for  Multiple  Delamination 

As  well  as  numerical  formulations,  we  have  continued  a  prior  line  of  research  (Andrews 
et  al.,  2004;  Andrews  et  al.,  2005;  Andrews  et  al.,  2006;  Cox  et  al.,  2004)  into  the 
applicability  of  advanced  beam  theory  methods.  This  effort  has  proven  a  good 
investment,  since,  in  many  cases  of  interest,  beam  theory  is  surprisingly  accurate  when 
appropriately  formulated  and  has  given  insight  into  deformations  and  energy  release  rates 
associated  with  multiple  delamination  cracking  that  may  not  have  been  nearly  as  clearly 
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highlighted  by  purely  numerical  work. 
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Figure  1.  (a)  X-ray  image  of  delamination  crack  growing  from  the  tip  of  a  sharp  notch 
(Spearing  and  Beaumont,  1992).  (b)  Prediction  in  which  delamination  and  splitting  cracks 
are  the  only  damage,  (c)  Prediction  in  which  delamination  and  splitting  cracks  are 
accompanied  by  distributed  softening  (due  to  microcracking)  in  the  laminate  plies. 


A  model  was  formulated  for  the  multiple  delamination  of  orthotropic  plates  subject  to 
static  and  dynamic  loading  and  deforming  in  cylindrical  bending  (Andrews  and  Massabo, 
2005).  The  model  was  based  on  a  particularization  of  Timoshenko  beam  theory  and 
accounts  for  cohesive  and  bridging  mechanisms  as  well  as  contact  and  friction  acting 
along  the  delamination  surfaces.  The  model  also  accounts  for  the  local  elastic 
deformation  at  the  delamination  tips,  which  give  rise  to  relative  rotations  between  the 
different  sub-beams  (root  rotations),  through  rotational  springs.  The  springs  arc  activated 
by  the  crack  tip  stress  resultants,  normal  and  shear  forces  and  bending  moments.  The 
compliances  of  the  springs  have  been  derived  numerically  for  a  wide  range  of  elastic 
constants.  Accounting  for  the  root  rotations  substantially  improves  the  accuracy  of  the 
one-dimensional  model,  leading  to  results  that  are  in  excellent  agreement  with  rigorous 
two-dimensional  solutions.  In  addition,  it  allows  for  closed  form  derivation  of  analytical 
expressions  for  the  energy  release  rate  and  stress  intensity  factors.  These  expressions  are 
an  extension  of  work  by  Suo  et  al.  and  Li  et  al.  to  account  for  shear,  near  tip  deformation 
and  multiple  delaminations.  They  depend  only  on  the  crack  tip  stress  resultants,  the 
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elastic  constants,  the  geometry  and  the  compliances  of  the  springs  describing  the  root 
rotations.  These  analytical  results  are,  we  feel,  an  important  contribution  to  the 
fundamental  development  of  beam  theory  for  fracture  studies. 

The  model  has  been  used  to  re-examine  single  delamination  specimens  of  standard 
fracture  tests.  The  validity  of  the  model  has  also  been  demonstrated  in  multiple 
delamination  systems,  where  mixed  mode  conditions  often  occur,  and  for  orthotropic 
materials  with  a  large  mismatch  in  elastic  parameters,  for  which  simplified  models 
neglecting  root  rotations  yield  inaccurate  predictions.  The  limitations  of  the  model,  as 
well  as  of  all  models  based  on  a  first  order  approximation  of  shear,  in  dealing  with  short 
process  zone  lengths  at  the  delamination  tip  where  cohesive  tractions  act  parallel  to  the 
delamination  surfaces  have  been  highlighted.  These  limitations  can  be  removed  using 
higher  order  theories. 

The  model  has  been  applied  to  investigate  the  interaction  effects  of  multiple 
delaminations  and  of  delaminations  and  regions  of  small  delamination  damage. 
Phenomena  of  amplification  and  shielding  of  the  energy  release  rate  and  modifications  in 
the  mode  ratio  have  been  observed.  The  different  regimes  of  behavior  have  been  defined 
for  a  system  of  two  delaminations  in  a  cantilever  beam. 

Finally,  in  a  phase  of  the  work  that  is  essential  for  the  current  program,  the  model  has 
been  used  to  investigate  the  effect  of  different  dynamic  loading  pulse  durations  on  the 
fracture  parameters  of  delaminated  plates.  Regimes  of  amplification  and  shielding  of  the 
energy  release  rate  due  to  dynamic  effects  have  been  identified  for  a  clamped  beam  with 
a  single  delamination  (Fig.  2)  and  the  results  have  been  summarized  using  shock  spectra. 
These  first  results  for  dynamic  delamination  will  be  used  in  our  future  work,  including 
numerical  simulations,  as  standards  for  limiting  cases. 


Figure  2.  Normalized  energy  release  rate  of  a  single  centrally  located  stationary  crack  in  a 
clamped  beam  specimen,  with  first  natural  period  7j  subject  to  different  pulse  durations,  Tp. 


The  question  of  whether  a  numerical  approach  or  beam  theory  is  most  appropriate  for  our 
further  studies  is  central  in  our  minds.  As  we  move  into  problems  with  more  complicated 
material  effects,  such  as  stitching  and  crack  face  friction,  as  well  as  specimens  or 
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components  with  non-trivial  geometry,  numerical  methods  will  obviously  be  favoured. 
However,  beam  theory  results  continue  to  prove  valuable  for  analyzing  limiting  cases,  for 
the  special  insight  they  can  yield  into  the  chief  factors  governing  structural  response. 

Dynamic  Fiber  Sliding 

In  a  collaboration  with  Professor  Ares  Rosakis  of  California  Institute  of  Technology,  we 
completed  a  detailed  study  of  the  physics  of  interfacial  sliding  in  the  prototypical 
problem  of  a  fiber  being  pulled  dynamically  out  of  a  matrix  (Yang  et  al.,  2006).  This 
work  was  initiated  under  a  prior  ARO  grant  (DAAD19-99-C-0042).  A  substantial 
amount  of  work  on  the  problem  and  the  writing  of  the  resulting  paper  were  then 
performed  without  funding,  after  that  grant  had  expired.  Under  the  new  grant  covered  by 
this  report  (W91 1NF-05-C-0073),  revisions  were  then  performed  in  August  and 
September,  2005,  in  response  to  some  interesting  ideas  raised  by  the  journal  reviewers. 
The  original  contributions  of  the  paper  include:  details  of  the  deformation  mechanics  of 
the  dynamic  pullout  problem;  confirmation  using  numerical  results  of  the  accuracy  of 
simple  analytical  shear  lag  solutions  (detailed  by  us  elsewhere  (Sridhar  et  al.,  2003));  and 
the  interesting  observation,  in  experiments  and  simulations,  of  indications  of  complex 
deformation,  perhaps  involving  chaotic  behaviour,  during  unloading.  The  contributions 
contained  specifically  in  the  revisions  made  under  this  contract  relate  to  the  inference  of 
information  about  interfacial  friction  from  our  experiments,  using  the  analysis  of  our 
numerical  simulations.  We  point  out  especially  that  the  disorder  found  upon  unloading  is 
insensitive  to  details  of  the  friction  law,  which  raises  an  essential  distinction  vis-a-vis  ill- 
posedness  and  instability  remarked  previously  in  the  literature,  which  depend  on  the 
friction  quite  sensitively.  Furthermore,  whereas  all  prior  studies  considered  time-uniform 
far-field  loading  rates,  our  observations  show  disorder  confined  to  the  unloading  regime  - 
the  loading  regime  is  smoothly  behaved.  The  possibility  that  unloading  creates 
qualitative  effects  in  deformation  fields  that  are  not  observed  for  rising  or  uniform 
loading  is  therefore  raised. 
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8.  Report  of  Inventions 


There  have  been  no  inventions  in  this  contract  to  date. 
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Abstract 

This  paper  deals  with  the  elastic  interaction  of  multiple  through-width  delaminations  in  laminated  plates  subject  to 
static  out  of  plane  loading  and  deforming  in  cylindrical  bending.  A  model  has  been  formulated  utilizing  the  classical 
theory  of  the  bending  of  beams  and  plates  and  accounting  for  non-frictional  contact  along  the  delamination  faces. 
Strong  interaction  effects  arise  between  the  delaminations  including  shielding  and  amplification  of  the  energy  release 
rate  and  modification  of  the  mode  ratio  as  compared  to  a  structure  with  only  a  single  delamination.  Such  behavior 
has  been  summarized  in  maps  that  completely  characterize  the  response  of  a  system  of  two  delaminations  in  a  cantilever 
beam.  The  quasi-static  propagation  of  the  system  of  delaminations  is  also  strongly  controlled  by  the  delamination  inter¬ 
actions,  which  lead  to  local  snap-back  and  snap-through  instabilities,  crack  arrest  and  crack  pull-along.  The  results 
show  similarity  to  those  for  cracked  infinite  bodies,  but  the  finite-thickness  of  the  plate  plays  an  important  role  and 
gives  rise  to  more  complex  behaviors.  The  stability  of  the  equality  of  length  of  a  system  of  n  delaminations  is  controlled 
by  their  spacing.  Finite  element  calculations  confirm  that  the  model  proposed  here  is  accurate,  except  when  the  differ¬ 
ence  in  the  length  of  the  interacting  delaminations  is  less  than  a  few  times  the  separation  of  their  planes. 

©  2005  Elsevier  Ltd.  All  rights  reserved. 

Keywords:  Multiple  delamination;  Shielding;  Amplification;  Crack  interaction;  Laminated  structures 


1.  Introduction 

Due  to  poor  interlaminar  properties,  laminated  fiber  reinforced  composites  are  susceptible  to  delamin¬ 
ations  caused  for  instance  by  manufacturing  errors,  edge  effects  or  by  such  in  service  situations  as  impact, 
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monotonic  and  cyclic  loading.  In  the  case  of  impact,  typically  many  delaminations  occur  between  layers  of 
different  fiber  orientation.  The  presence  of  the  delaminations  is  often  undetectable  on  the  surface  and  may 
significantly  reduce  the  stiffness  and  load  carrying  capacity  of  the  structure  (Pavier  and  Clarke,  1995).  When 
a  load  is  applied  that  causes  sufficient  interlaminar  stresses,  the  delaminations  may  grow,  often  catastroph¬ 
ically,  separating  the  structure  into  two  or  more  pieces  and  causing  further  reduction  in  stiffness  and/or  fail¬ 
ure.  While  the  behavior  of  structures  in  the  presence  of  a  single  delamination  has  been  widely  studied  since 
the  early  work  of  Kanninen  (1973,  1974)  and  others  (see  Massabo  and  Cox,  1999;  Sridhar  et  al.,  2002  for 
work  of  the  authors),  the  problem  of  the  multiple  delamination  of  laminated  structures  is  not  yet  fully 
understood.  This  paper  examines  the  interaction  of  multiple  delaminations  in  plates  subject  to  static  out 
of  plane  loading  and  the  effects  this  interaction  has  on  the  fracture  behavior  and  structural  response  of 
the  laminate,  using  a  cantilever  plate  as  a  case  study. 

An  important  interaction  effect  is  contact  that  may  occur  between  the  delaminated  plies.  When  there  are 
multiple  delaminations,  studies,  such  as  those  dealing  with  structures  subject  to  in  plane  loadings  (Larsson, 
1991;  Suemasu,  1993),  have  shown  that  extensive  contact  may  occur  along  the  delaminations  faces.  Contact 
significantly  affects  the  critical  energy  release  rate  for  the  extension  of  a  delamination.  Contact  also  intro¬ 
duces  regions  in  which  friction  may  be  important.  The  presence  of  contact  will  be  shown  later  not  to  be 
limited  to  delaminated  structures  subject  to  in-plane  loading. 

The  effect  of  the  interaction  of  the  delaminations  on  the  energy  release  rates  is  another  important  phe¬ 
nomenon.  Larsson  (1991),  using  a  model  based  on  the  theory  of  bending  of  plates  to  study  delaminated 
plates  under  in-plane  loading,  observed  a  discontinuity  in  the  energy  release  rate  of  a  delamination  when 
its  length  is  equal  to  the  length  of  other  delaminations  of  the  system.  This  discontinuity  was  thought  to 
be  an  anomalous  product  of  the  assumptions  of  plate  theory;  however,  this  paper  will  show  that  the  dis¬ 
continuity  well  approximates  the  actual  behavior.  Zheng  and  Sun  (1998)  showed  that  the  effect  of  the  inter¬ 
actions  between  delaminations  on  the  energy  release  rates  depends  on  the  through-thickness  distance 
between  the  delaminations.  For  the  structure  and  crack  configurations  that  they  studied,  they  also  observed 
that  the  presence  of  a  smaller  delamination  has  little  effect  on  the  energy  release  rate  of  a  delamination, 
while  the  presence  of  a  longer  delamination  may  induce  strong  amplification  or  shielding.  This  paper  will 
show  that  in  more  general  crack  configurations,  the  presence  of  short  delaminations  may  also  strongly  af¬ 
fect  the  response. 

The  interaction  effects  also  influence  the  stability  of  the  equality  of  length  of  a  system  of  n  delaminations. 
As  it  will  be  shown  later,  the  delaminations  of  a  stable  system  will  tend  to  grow  together,  leading  to  an 
increased  capacity  of  the  system  to  absorb  energy  and  a  more  ductile  structural  response.  An  unstable  sys¬ 
tem  of  delaminations  will  have  more  localized  delamination  growth  of  only  one  or  a  small  number  of  del¬ 
aminations,  leading  to  reduced  energy  absorption  and  a  less  ductile  response.  Suemasu  and  Majima  (1996) 
showed  that  an  axisymmetric  system  of  equally  spaced,  equal  size,  penny-shaped  delaminations  in  a 
clamped  circular  plate  subject  to  a  concentrated  static  force  is  stable.  More  general  conclusions  on  the  sta¬ 
bility  of  systems  of  delaminations  with  unequally  spaced  cracks  will  be  presented  in  the  following. 

The  model  proposed  in  this  paper  is  based  on  the  theory  of  bending  of  beams  and  plates  and  examines  a 
multiply  delaminated  plate  subject  to  static,  out  of  plane  loading  and  deforming  in  cylindrical  bending. 
Non-frictional  contact  along  the  delamination  surfaces  is  accounted  for  utilizing  three  approximations, 
two  of  which  allow  closed  form  solutions  of  the  problem.  The  energy  release  rate  and  the  stress  intensity 
factors  at  each  delamination  crack  tip  are  determined  by  the  application  of  methods  developed  for  beams 
and  plates  with  a  single  delamination.  The  results  are  quantified  for  the  case  of  a  cantilever  beam  with  mul¬ 
tiple  edge  delaminations.  Shielding  and  amplification  of  the  energy  release  rates  of  the  cracks  and  the  inter¬ 
action  effects  on  the  macro  structural  response  of  the  plate  are  examined.  The  model  is  validated  through 
finite  element  analyses. 
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2.  Theoretical  model 

2.1.  Model  assumptions 

The  problem  under  consideration  is  a  thin  laminated  plate  of  thickness  h  with  n  through  width  delam¬ 
inations.  A  system  of  Cartesian  coordinates  x-y-z  is  introduced.  Each  delamination,  with  index  i,  is  arbi¬ 
trarily  distributed  across  the  thickness,  with  through-thickness  position  yf  and  length  ak  The  plate  is  subject 
to  out  of  plane  loading.  An  exemplary  plate,  with  an  applied  surface  pressure  q ,  is  shown  in  Fig.  la. 

The  material  comprising  the  plate  is  assumed  to  be  homogeneous,  isotropic  and  linear  elastic,  with  elas¬ 
tic  modulus  E  and  Poisson’s  ratio  v.  The  model  is  applicable  to  all  laminates  composed  of  isotropic  layers  of 
the  same  material.  With  modifications  to  account  for  different  material  properties  in  the  through  thickness 
direction,  the  model  can  also  be  applied  to  quasi-isotropic  laminates  and  specially  orthotropic  laminates 
(i.e.  laminates  composed  of  layers  whose  principal  material  directions  are  aligned  with  the  reference  axes) 
with  a  large  number  of  layers. 

Only  small  deformations  of  the  plate  are  considered,  and  plane  strain  conditions  are  assumed  parallel  to 
the  x-y  plane  so  that  the  plate  deforms  in  cylindrical  bending.  Therefore,  the  governing  equations  of  the 
plate  correspond  to  those  of  a  beam  with  a  reduced  Young’s  modulus  E  =  E/(  1  —  v2),  where  E  is  the  lon¬ 
gitudinal  modulus.  In  the  following  analysis,  the  designation  of  beam  and  plate  are  identical.  The  classical 
theory  of  bending  of  beams,  which  neglects  shear  deformation,  is  utilized  to  determine  the  response  of  the 
delaminated  structure. 

The  plate  shown  in  Fig.  la,  is  decomposed  into  multiple  uncracked  beam  segments,  defined  by  the  coor¬ 
dinates  of  the  crack  tips  and  continuity  conditions  are  applied  at  the  cross  sections  separating  beam  seg¬ 
ments.  For  purposes  of  notation,  the  beam  segments  are  numbered  from  top  to  bottom  and  left  to  right 
in  the  structure,  and  delaminations  are  numbered  from  top  to  bottom.  The  height  of  beam  segment  k  is 
hk  and  its  cross  sectional  moment  of  inertia  and  area  are  Ik  and  Ak ,  respectively.  The  generalized  displace¬ 
ments  of  the  centroidal  axis  of  beam  segment  k  are  the  axial  and  transverse  displacements  uk  and  wk,  and 
the  bending  rotation  cpk.  The  stress  resultants  per  unit  width  are  axial  force  Nk,  shear  force  Vk  and  bending 
moment  Mk  (Fig.  lb). 
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Fig.  1.  (a)  Exemplary  plate  with  multiple  through  width  cracks  subject  to  transverse  loading,  (b)  Equilibrium  of  a  beam  segment, 
showing  stress  resultants  and  contact  pressures,  (c)  Crack  tip  element  separating  beam  segments  identified  by  j,  k  and  k+  1.  The 
dashed  line  is  the  path  used  for  calculating  the  J  integral. 
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Assuming  no  shear  deformation,  yk  =  0,  the  compatibility  equations  for  the  beam  segment  k  are 

=  £k  =  UktX,  Kk  =  (pkj,  (1) 

where  ek  is  the  axial  deformation  and  Kk  is  the  bending  curvature  and,  x  denotes  differentiation  with  respect 
to  the  longitudinal  coordinate  x.  The  constitutive  equations  are 

EIkKk  El  kWkxx  •>  Nk  EAkSk  EA.kUkx^  (2) 

Equilibrium  of  the  beam  segment  is  given  by: 

Mk,x  —  Vk  0,  EkjX  +  Pk,k+ 1  Pk-lji  Nk,x  0}  (3) 

where  pk,k+\  and  Pk-\,k  are  either  externally  applied  pressures  or  contact  pressures  acting  on  the  lower  and 
upper  surfaces  of  the  segment  (Fig.  lb). 

At  the  zth  crack  tip,  Fig.  lc,  with  coordinate  x  =  x,  separating  the  beam  segments j,  k  and  k+  1,  the  kine¬ 
matic  and  static  continuity  conditions  are 

Mj  T  ~  (hj  hk)wjjX  uk:  Uj  —  (hj  hk+\)wjjX  uk+\, 

Wj  =Wk  =  Wk+ 1, 

WPC  =  =  Wk+1„ 

Nj=Nk  +  Nk+u  [  j 

Mj  =  Mk  +  Mk+ 1  —  -  (hj  —  hk)Nk  +  -  (hj  —  hk+i)Nk+i, 

Vj  =  Vk  +  Vk+1. 

The  bending  rotations  of  the  beam  segments  at  the  crack  tips,  often  referred  to  as  root  rotations,  have  been 
assumed  in  (4)  to  be  identical  to  one  another. 

In  the  plate  shown  in  Fig.  la,  depending  on  the  loading  conditions  and  the  geometry  of  the  delamination 
system,  contact  along  the  delamination  surfaces  may  occur.  This  contact  is  assumed  to  be  non-frictional, 
which  allows  free  sliding  along  the  surfaces  of  a  delamination,  and  is  represented  by  the  contact  pressures, 
pk,k+  i  and  Pk-\,k •  The  presence  of  contact  is  simulated  in  three  ways. 

The  first  method  is  to  assume  that  the  deflections  of  the  beams  in  the  delaminated  regions  are  the  same. 
Interpenetration  along  the  delamination  surfaces  is  then  avoided;  however,  the  beam  segments  in  the  del¬ 
aminated  region  are  also  constrained  from  separating  from  each  other,  thus  preventing  any  opening  along 
the  delamination  surfaces  that  may  occur.  This  is  referred  to  as  the  constrained- contact  model  and  intro¬ 
duces  the  constraint  equation  wk  =  wk+\  between  the  delamination  surfaces  of  beam  segments  k  and  k+  1. 

The  second  method  is  to  assume  no  contact  interaction  between  the  beam  segments  in  the  delaminated 
region.  This  allows  opening  along  the  delamination  surfaces;  however,  interpenetration  can  occur, 
wk  ^  wk+ 1.  This  is  referred  to  as  the  unconstrained- contact  model. 

The  final  method  is  to  consider  elastic  contact  along  the  delamination  surfaces.  The  elastic  contact  is 
approximated  using  a  Winkler  foundation  of  linear  springs,  which  represent  the  through-thickness  stiffness 
of  the  contacting  beam  segments  and  act  to  resist  interpenetration  along  the  delamination  surfaces.  The 
stiffness  of  the  springs  for  a  contacting  beam  segment  k  is  determined  by  considering  one  dimensional  ten¬ 
sion  or  compression  of  a  beam  element  of  height  hj 2.  This  model  neglects  shear  deformations  in  the  beam 
segment,  which  would  lead  to  coupling  of  the  springs  (see  Kerr  (1964)  for  an  in-depth  discussion  of  foun¬ 
dation  models),  in  accordance  with  what  has  been  assumed  initially.  The  stiffness  of  the  layer  of  springs 
representing  contact  between  two  beam  segments,  k  and  k  +  1,  is  given  by: 

2  Et 


h,k+ lM  =  H(wk(x )  -  wk+x(x)) 


hk  +  hk+ 1 


(5) 
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where  ET  is  the  through-thickness  modulus  of  the  material  (ET  =  E  for  an  isotropic  material)  and  Hf)  is  the 
Heaviside  step  function,  //(£)  =  { 1,£  >  0;0,£  <  0},  which  ensures  that  the  springs  resist  only  interpenetra¬ 
tion  and  do  not  resist  opening  between  beam  segments.  The  pressure  exerted  at  the  bottom  and  top  of  beam 
segments  k  and  k  +  1  then  follows  as: 

Pk,k+ 1 M  =  kk,k+ 1  (x)(wk(x)  -  wk+i  (x)).  (6) 

This  method  is  referred  to  as  the  spring-contact  model.  Since  the  regions  of  contact  and  opening  are  gen¬ 
erally  not  known  a  priori,  the  problem  in  this  case  becomes  nonlinear. 

It  will  be  shown  in  the  following  that  the  constrained-  and  unconstrained- contact  models  define  upper  and 
lower  bound  solutions  of  the  spring-contact  model.  Both  models  can  be  recovered  from  the  spring-contact 
model  by  taking  the  appropriate  limit  of  the  spring  stiffness  in  Eq.  (5).  For  the  constrained-contact  model, 
the  limit  kk^+ 1  — ►  oo  results  in  the  constraint  equation  already  defined,  wk  =  wk+ 1.  For  the  unconstrained- 
contact  model,  the  limit  as  kk  k+\  =  0  results  in  pk,k+\{x)  =  0,  namely  the  absence  of  any  contact  pressure. 
The  constrained- and  unconstrained-contact  models,  while  less  accurate  than  the  spring- contact  model,  have 
the  advantage  of  leading  to  closed  form  solutions  of  the  problem  in  cases  that  require  a  numerical  solution 
when  treated  with  the  spring-contact  model.  The  solutions  obtained  with  the  simplified  models  are  qualita¬ 
tively  similar  to  those  of  the  spring-contact  model  and  they  can  lead  to  insightful  conclusions. 

The  proposed  model  relies  on  the  three  simplifying  assumptions  of  zero  shear  deformation,  absence  of 
relative  beam  root  rotations  at  the  crack  tips  and  non-frictional  contact.  The  influence  of  shear  deformation 
on  the  solution  of  single  delamination  problems  is  known  to  affect  only  quantitative  details  of  the  solution 
and  to  be  negligible  if  the  delamination  is  sufficiently  long.  In  multiple  delamination  problems  a  stronger 
effect  is  expected  when  the  delamination  tips  are  close.  However,  current  studies  (Andrews  et  al.,  2005) 
show  that  this  effect  is  not  strong  in  the  absence  of  contact  near  the  delamination  tips  and  is  negligible  com¬ 
pared  to  the  effect  of  unequal  beam  rotations  at  the  delamination  tips  in  the  presence  of  contact. 

The  influence  of  the  assumption  of  equal  beam  rotations  at  the  crack  tips  has  been  shown  previously  to 
affect  only  quantitative  details  of  the  solution  for  single  delaminations,  leading  to  limited  underestimation 
of  the  compliance  of  the  system  and  the  energy  release  rate  for  mode  I  fracture  problems  and  to  negligible 
effects  in  mode  II  fracture  problems.  In  the  presence  of  multiple  delaminations  the  effect  of  the  assumption 
is  expected  to  be  stronger  due  to  the  presence  of  contact  between  the  crack  faces.  Different  methods  for 
correcting  root  rotations  have  been  proposed  since  the  early  work  of  Kanninen  (1973,  1974)  (see  for  in¬ 
stance  Williams,  1989;  Sun  and  Pandey,  1994,  Pandey  and  Sun,  1996;  and  Wang  and  Qiao,  2004).  All  meth¬ 
ods,  however,  complicate  the  solution  of  the  problem  and  as  the  crack  interaction  effects  examined  in  this 
work  will  arise  independently  of  the  quantitative  details  of  crack  face  phenomena,  a  simpler  level  of  approx¬ 
imation  has  been  utilized  for  this  work. 

Accounting  for  the  presence  of  friction  between  the  crack  surfaces  does  not  complicate  the  model  greatly 
and  could  lead  to  interesting  alterations  of  the  results  presented  in  this  paper  for  certain  geometries.  This 
effect  is  not  studied  here  and  is  the  topic  of  current  work  (Andrews  et  al.,  2005).  The  assumptions  of  the 
model  and  their  implications  will  be  discussed  further  in  Section  5. 

2.2.  Model  solution 

An  arbitrary  section  of  the  beam  shown  in  Fig.  la  may  be  intersected  by  m  of  the  n  delaminations  in  the 
beam.  Using  the  compatibility,  equilibrium  and  constitutive  equations  defined  above,  the  governing  differ¬ 
ential  equations  for  the  beam  segments  intersected  by  the  section  are 

EIkwktXXXX  +  pkk+l  =  pk_lk,  k=l,...,m+l, 

EAkuktXX  =  0,  k  =  1, . . .  ,m  +  1, 


(7) 

(8) 
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where  the  p's  are  contact  pressures  between  the  beams,  except  for  p0  i  and  pm+ i?m+2  which  are  externally 
applied  pressures  on  the  external  surfaces  of  the  section.  Eq.  (8),  referring  to  the  axial  displacements,  are 
not  coupled  and  general  solutions  can  be  derived  for  all  k's.  Similarly,  Eq.  (7),  referring  to  the  transverse 
displacements,  are  not  coupled  and  can  be  solved  analytically  within  the  approximation  of  unconstrained- 
contact.  Full  coupling  derives  from  the  constraint  equations,  =  w^+i,  of  the  constrained-contact  model 
that  again  allows  for  a  closed  form  general  solution. 

For  the  spring- contact  model,  on  the  other  hand,  the  p's  are  given  by  Eq.  (6),  which  results  in  a  system 
of  m  +  1  coupled  differential  equations,  Eq.  (7).  Because  the  springs  representing  contact  do  not  act  in 
tension,  contact  may  not  occur  between  all  of  the  beam  segments  intersected  by  this  section.  The  presence 
of  traction-free  crack  surfaces  simplifies  the  system  of  coupled  Eq.  (7)  so  that  a  set  of  subsystems  char¬ 
acterized  by  a  lower  number  of  coupled  variables  and  equations  can  be  defined.  The  limit  cases  are  those 
of  full  coupling  of  the  wk  for  k  — 1, . .  .,ra  +  1,  when  contact  exists  between  all  segments,  and  no  coupling 
when  all  crack  surfaces  are  opened.  The  coupled  differential  equations  are  linear  with  constant  coeffi¬ 
cients  and  the  characteristic  algebraic  equations  for  the  sub-systems  can  always  been  found.  The  general 
closed  form  solution  of  a  subsystem  composed  of  one  free  segment  and  of  2  or  3  contacting  beam  seg¬ 
ments  is  shown  in  Appendix  A.  A  subsystem  of  more  than  three  beam  segments  generally  requires 
numerical  solution.  General  solutions  of  the  subsystems  for  each  section  of  the  beam  in  all  possible  states 
of  opening  and  contact  can  then  be  determined  and  used  in  the  iterative  procedure  to  define  the  regions 
of  contact  and  opening. 

The  solution  for  the  whole  structure  is  determined  by  applying  the  kinematic  and  static  boundary  and 
continuity  conditions  to  define  the  constants  of  integration,  Q,  of  the  general  solutions  of  systems  (7)  and 
(8).  If  the  unconstrained-contact  or  the  constrained-contact  model  is  used,  the  problem  can  then  be  solved  in 
closed  form.  If  the  unconstrained-contact  model  predicts  no  interpenetration  between  the  crack  surfaces,  the 
solution  is  exact.  In  addition,  the  regions  of  crack  face  interpenetration  predicted  through  the  uncon¬ 
strained-contact  model  give  a  first  approximation  of  the  regions  where  contact  is  expected  to  occur  within 
the  spring-contact  model. 

If  the  spring-contact  model  is  used,  the  regions  where  crack  face  contact  and  opening  take  place  are  un¬ 
known  a  priori  and  must  be  determined  by  an  iterative  solution  process.  The  iterative  process  is  initiated  by 
assuming  an  initial  state  of  contact  along  the  delamination  surfaces.  The  beam  segments  formed  by  the  divi¬ 
sions  at  the  crack  tips  are  further  subdivided  at  all  coordinates  where  there  is  a  change  in  contact  state  along 
a  delamination  surface.  At  these  coordinates,  additional  continuity  conditions  are  imposed  between  subdi¬ 
visions  j  and  k ,  which  correspond  to  Eqs.  (4)  with  the  contributions  from  beam  segment  k  +  1  removed.  A 
system  of  algebraic  equations,  corresponding  to  the  boundary  and  continuity  conditions  at  all  divisions  of 
the  beam  and  the  general  solutions  of  the  governing  Eqs.  (7)  and  (8),  is  constructed  based  on  this  assumed 
state  of  contact  and  is  solved  numerically  for  the  unknown  constants  of  integration,  Q  Utilizing  the  dis¬ 
placement  solution,  updated  regions  of  delamination  surface  interpenetration  and  opening  are  determined 
for  all  delaminations  in  the  system.  The  updated  regions  of  interpenetration  are  then  assumed  to  be  in  con¬ 
tact  and  a  new  system  of  algebraic  equations  is  constructed  and  solved  for  the  new  C/s.  A  grid  of  points  is 
introduced  along  the  axis  of  each  beam  segment  and  the  convergence  of  the  solution  is  checked  by  consid¬ 
ering  the  norm  of  the  transverse  displacements  at  these  locations.  The  process  is  repeated  until  convergence 
of  the  solution  to  a  specified  tolerance.  For  all  cases  examined  in  this  paper,  convergence  of  the  solution  has 
been  reached  easily  with  only  a  limited  number  of  iterations. 

2.3.  Energy  release  rate  and  mode  decomposition 

It  is  assumed  that  the  n  delaminations  in  the  structure  are  constrained  to  propagate  along  the  low  tough¬ 
ness  fracture  paths  defined  by  their  planes.  The  energy  release  rate  for  the  individual  coplanar  extension  of  a 
delamination  tip  in  this  system  is  determined  by  application  of  the  /-Integral  (Rice,  1968)  along  a 
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path  around  the  delamination  tip.  The  expression  for  the  /-integral  for  crack  tip  i,  at  the  coordinate  xh  sep¬ 
arating  beam  segments  j,  k  and  k+  1  (Fig.  lc),  is 


1 

(m\  N\  M2k+l 

M) 

N2j\ 

2  1 

\Eh  EAk  EIk+1 

EAk+x 

~eTj 

£ 

1 

This  expression  is  identical  to  the  expression  valid  for  beams  with  a  single  delamination,  as  is  explained  in  the 
following.  Each  crack  tip  in  the  beam  shown  in  Fig.  la  can  be  extracted  from  the  structure  as  an  equivalent 
beam  with  a  single  delamination  whose  upper  and  lower  surfaces  are  defined  either  by  the  surfaces  of  the 
plate  or  by  other  delaminations,  Fig.  lc.  This  equivalent  single  beam  is  subject  to  end  forces  and  moments 
as  well  as  possible  contact  pressures  on  its  upper  and  lower  surfaces.  The  path  for  the  /-integral,  shown  in 
Fig.  lc,  is  taken  at  the  crack  tip  cross  sections  immediately  preceding  and  following  the  crack  tip  and  along 
the  upper  and  lower  surfaces  of  the  equivalent  beam.  For  this  path,  the  contact  pressures  on  the  upper  and 
lower  surfaces  do  not  enter  the  expression  for  the  /-integral.  If  there  is  another  delamination  tip  at  the  same 
coordinate  xt  the  energy  release  rate  is  determined  as  the  limit  for  Axt  tending  to  zero  of  Eq.  (9a)  applied  to  a 
system  in  which  the  length  of  the  delamination  of  interest  has  been  increased  (the  sign  of  Ax,-  is  chosen  cor¬ 
responding  to  an  increase  in  delamination  length)  by  incrementing  the  position  of  the  crack  tip  to  x,-  +  Ax,. 
Expression  (9a)  is  only  valid  when  the  rotations  of  the  beams  at  the  crack  tips  are  assumed  to  be  the  same, 
and  must  be  modified  accordingly  if  this  assumption  is  relaxed  (see  Section  5).  Similarly,  the  equation  could 
be  modified  to  include  the  contributions  due  to  the  shear  deformations  along  the  beams. 

The  above  approach  cannot  be  applied  to  analyze  the  simultaneous  propagation  of  m  delamination  tips 
that  have  the  same  coordinate  xm.  Instead,  the  energy  release  rate  in  this  case  is  conveniently  determined  by 
the  variation  of  the  total  potential  energy  with  respect  to  simultaneous  unit  extension  da.  The  energy  release 
rate  for  each  of  the  delaminations  under  this  condition  is 


1  d  W 
m  da 


(9b) 


where  W  is  the  total  potential  energy  of  the  system.  The  value  of  da  should  be  chosen  such  that  it  corre¬ 
sponds  to  an  increase  in  length  of  the  delaminations.  For  the  coordinate  system  described  in  Fig.  1,  the 
extension  of  the  left  tips  of  delaminations  corresponds  to  da  =  —  dxm  and  the  extension  of  the  right  tips 
of  the  delaminations  corresponds  to  da  =  dxm.  The  variation  of  the  total  potential  energy  is  also  used  to 
analyze  the  simultaneous  propagation  of  m  delamination  tips  that  do  not  share  the  same  coordinate. 

The  conditions  at  the  crack  tips  in  the  general  beam  of  Fig.  la  are  generally  mixed  mode.  The  energy  re¬ 
lease  rate  calculated  using  Eqs.  (9a)  and  (9b)  includes  both  the  mode  I  and  mode  II  components,  ^  and  ^n,-. 
Separation  of  the  modes  of  fracture  is  accomplished  by  using  the  method  proposed  by  Suo  (1990),  Suo  and 
Hutchinson  (1990)  and  Hutchinson  and  Suo  (1992).  They  derived  analytical  expressions  for  the  mixed  mode 
stress  intensity  factors  for  a  beam  with  a  single  crack,  with  total  energy  release  rate  given  by  Eq.  (9a),  in  terms 
of  the  stress  resultants,  bending  moments  and  normal  forces,  at  the  cross  sections  immediately  preceding  and 
following  the  crack  tip,  and  geometrical  parameters.  These  expressions  depend  on  an  additional  parameter 
that  can  be  derived  with  a  rigorous  two-dimensional  solution  of  the  problem.  In  the  model  proposed  in  this 
paper,  the  method  is  applied  to  the  equivalent  beam  with  a  single  crack  shown  in  Fig.  lc  with  energy  release 
rate  given  by  Eq.  (9a)  and  is  thus  valid  for  each  crack  tip  in  the  system. 


3.  Cantilever  beam  with  n  equal  length  delaminations 

The  first  study  problem  is  a  cantilever  beam  of  length  L  with  n  through  width,  arbitrarily  spaced  edge 
delaminations  (Fig.  2)  subjected  to  a  static  concentrated  out  of  plane  force  P  at  its  free  end.  This  simple 
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Fig.  2.  Cantilever  beam  with  n  edge  cracks  subject  to  a  transverse  load  P. 


structure  provides  an  analytical/semi-analytical  description  of  how  multiple  delaminations  interact.  This 
problem  yields  results  that  are  identical  to  those  for  an  edge  cracked  simply  supported  beam,  loaded  by 
an  out  of  plane  force  at  its  mid  span.  Solutions  for  simply  supported  or  fixed  plates  with  mid-span,  out 
of  plane  loading  can  be  simply  obtained  by  modifying  the  boundary  conditions  at  x  =  0  and  x  —  L.  The 
solution  can  also  be  easily  extended  to  circular  axisymmetric  plates  loaded  by  concentrated  out  of  plane 
forces  at  the  centers. 

In  this  problem  no  axial  forces  are  developed;  thus  the  equations  governing  the  axial  displacements  are 
ignored.  The  boundary  conditions  for  this  problem  are 

wo  =  0,  wo)JC  =  0  at  x  =  0, 

j+n 

P=EVk'  Mk  =  0,  k  —  j, . . .  ,j  +  n,  (10) 

k=j 

Wk  =  w*+i,  k  =  j1...J  +  n-  1  at  x  =  Z, 

where  j  is  the  index  of  the  uppermost  beam  segment  at  the  free  end  of  the  beam;  and  as  before,  the  sub¬ 
scripts  define  the  number  of  the  beam  segments.  Contact  of  the  beam  segments  under  the  load  point  is 
approximated  by  assuming  that  their  deflections  at  the  free  end  are  equal.  This  assumption  is  exact  for 
the  unconstrained- contact  model  and  accurate  for  the  spring-contact  model  if  the  delaminations  are  long. 
The  static  and  kinematic  continuity  conditions  at  each  crack  tip  are  given  by  Eq.  (4). 

In  the  simplest  configuration,  the  delaminations  in  the  beam  shown  in  Fig.  2  are  of  the  same  length  a.  In 
this  configuration,  the  beam  is  divided  up  into  n  +  2  beam  segments,  n  +  1  segments  forming  the  cracked 
region,  and  one  segment  forming  the  uncracked  region.  Each  crack  tip  is  located  at  the  same  coordinate 
x  =  L  —  a.  The  continuity  conditions  (4)  are  modified  as  follows: 

w0  =  wk,  W0)X  =  wkjX,  k  =  1, . . . ,  n  +  1, 

w+l  n+ 1  n\) 

m0  =  J2m^  vo  =  J2Vk- 

k=l  k=  1 

Combining  the  general  solutions  for  each  beam  segment  for  the  three  contact  models  allows  this  system  of 
A(n  +  2)  boundary  and  continuity  equations  to  be  solved  in  closed  form  (Appendix  B).  The  three  methods 
yield  identical  solutions,  which  shows  that  the  deflections  of  the  beam  segments  in  the  cracked  region  of  the 
beam  are  identical  and  there  is  no  contact  or  opening  along  the  crack  faces. 

3.1.  Energy  release  rates  for  simultaneous  propagation  of  equal  length  delaminations 

The  energy  release  rate  for  each  delamination  when  the  delaminations  are  assumed  to  propagate  simul¬ 
taneously  is  conveniently  determined  by  the  variation  of  the  total  potential  energy  W  due  to  a  unit  prop¬ 
agation  of  the  cracks,  Eq.  (9b): 
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1  dW  _  1  P2a2  ( 
n  d a  2 n  EIq  l  0 


(12) 


where  /0  is  the  moment  of  inertia  of  the  intact  portion  of  the  beam  and  /&  is  the  moment  of  inertia  of  the 
beam  segment  in  the  cracked  portion  of  the  beam.  The  next  sections  will  show  that  for  n  equally  spaced 
delaminations  the  energy  release  rate  for  simultaneous  propagation  is  higher  than  the  energy  release  rate 
corresponding  to  the  propagation  of  one  of  the  delaminations  of  the  system.  On  the  other  hand,  when 
the  delaminations  are  unequally  spaced,  the  response  is  controlled  by  the  transverse  position  of  the 
delaminations. 

Eq.  (12)  shows  that  if  the  delaminations  are  uniformly  spaced  across  the  entire  cross  section,  then  as  the 
number  of  delaminations  increases,  the  energy  release  rate  per  delamination  also  increases.  If  on  the  other 
hand,  the  cracks  are  confined  to  a  band,  then  as  the  number  of  cracks  within  the  band  increases  the  energy 
release  rate  per  crack  decreases.  The  latter  case  could  represent  damage  localization  in  the  beam.  If  a  crack 
growth  criterion  based  on  the  total  energy  release  rate  is  applied,  the  first  case  would  correspond  to  a  de¬ 
crease  of  the  critical  load  on  increasing  the  number  of  cracks,  and  the  second  case  to  an  increase  of  the 
critical  load.  However,  for  the  second  case,  since  the  cracks  are  not  equally  spaced,  the  assumption  of 
simultaneous  propagation  may  not  be  valid  and  the  critical  load  could  be  defined  by  the  propagation  of 
only  one  of  the  cracks  of  the  system. 


3.2.  Stability  of  the  equality  of  length  in  a  system  of  equal  length  delaminations 

The  equality  of  length  will  be  said  to  be  stable  in  a  system  of  equal  length  delaminations  if  the  system 
recovers  the  condition  of  equal  length  after  the  lengths  of  one  or  more  cracks  are  perturbed.1  It  will  be 
shown  later  that,  for  a  fixed  number  of  delaminations,  the  energy  absorption  is  higher  if  the  equality  of 
length  is  stable  than  if  it  is  not. 

3.2.1.  Equally  spaced  delaminations 

This  special  case  was  previously  examined  by  Suemasu  and  Majima  (1996)  who  investigated  the  stability 
of  the  equality  of  length  of  an  axisymmetric  system  of  equally  spaced,  penny-shaped  delaminations  in  a 
clamped  circular  plate  subject  to  a  concentrated  out-of-plane  load.  The  study  was  performed  under  the  sim¬ 
plifying  kinematic  assumption  of  constrained-contact.  This  assumption  is  relaxed  here,  and  the  stability  of 
the  equality  of  length  is  studied  using  the  spring-contact  model.  The  problem  is  greatly  simplified  by  noting 
that  in  the  limit  of  the  delaminations  having  the  same  length,  the  solution  of  the  spring- contact  model  is 
identical  to  that  of  the  unconstrained-contact  model.  Closed  form  solutions  can  therefore  be  found,  which 
are  detailed  in  Appendix  C  and  summarized  below. 

An  ideal  system  of  equal  length,  equally  spaced  delaminations  propagates  when  the  energy  release  rate, 
given  by  Eq.  (12),  equals  the  fracture  energy  of  the  material.  If  a  positive  perturbation  of  the  length  of  crack 
i  in  the  system  is  considered,  then  the  energy  release  rate  for  the  propagation  of  that  crack  is  lower  than  that 
corresponding  to  the  simultaneous  propagation  of  the  remaining  cracks.  Thus,  using  a  fracture  criterion 
based  on  the  total  energy  release  rate,  the  remaining  delaminations  will  grow  up  to  the  length  of  the  per¬ 
turbed  delamination  and  equality  of  length  will  be  restored.  Conversely,  if  a  negative  perturbation  of  the 
length  of  crack  i  is  considered,  then  the  energy  release  rate  of  crack  i  is  higher  that  of  the  remaining  cracks: 
the  system  is  again  stable.  Since  a  negative  perturbation  of  one  crack  is  identical  to  a  positive  perturbation 


1  The  term  stability  is  used  here  to  refer  to  changes  in  the  crack  geometry  and  not  the  question  of  whether  the  crack  growth  is  quasi¬ 
static  (stable)  or  dynamic  (unstable). 


864 


M.  G.  Andrews  et  al.  /  International  Journal  of  Solids  and  Structures  43  ( 2006 )  855-886 


of  the  remaining  cracks  and  vice  versa,  conclusions  on  the  effect  of  positive  or  negative  perturbations  of  one 
crack  will  also  hold  for  perturbations  of  a  generic  number  of  cracks. 

The  same  conclusions  are  reached  if  the  simplifying  kinematic  assumption  of  constrained  contact  is 
used.  The  analysis,  which  is  similar  to  that  presented  in  the  Appendix  C,  is  not  presented  here.  Suemasu 
and  Majima  (1996)  obtained  the  same  results  for  axisymmetric  clamped  circular  plates. 

3.2.2 .  Unequally  spaced  delaminations 

A  system  of  unequally  spaced,  equal  length  delaminations  does  not  always  grow  self-similarly,  even  in 
the  absence  of  length  perturbations.  The  question  of  stability  of  the  equality  of  length  is  controlled  by 
the  through-thickness  position  of  the  delaminations.  The  map  of  Fig.  3  refers  to  a  two-crack  system  and 
has  been  constructed  by  using  the  spring- contact  model.  Similar  results  are  obtained  using  the  uncon¬ 
strained-contact  model.  Fig.  3  shows  regions  in  which  the  energy  release  rates  for  the  propagation  of  one 
of  the  two  cracks,  and  ^L,  and  for  their  simultaneous  propagation,  ^B,  have  different  ordering.  The 
map  is  a  function  of  the  through-thickness  positions  of  the  two  cracks.  The  upper  crack  is  located  at  a  dis¬ 
tance  h 3  from  the  upper  surface  of  the  beam  and  the  lower  crack  is  located  at  a  distance  h5  from  the  lower 
surface  of  the  beam,  as  shown  in  the  inset  in  the  figure. 

The  map  shows  three  regions:  region  I  in  which  >  max(^L,^B);  region  II  in  which 
^B  >  max(^u,  ^L);  and  region  III  in  which  >  max(^U5  ^B).  Thus  in  regions  I  and  III,  only  one  crack 
will  propagate  when  critical  conditions  are  satisfied,  while  in  region  II,  both  cracks  propagate  together. 
Furthermore,  in  regions  I  and  III,  the  system  is  more  unstable  (increasing  difference  between  the  highest 
and  next  highest  values  of  and  ^B)  under  positive  perturbations  of  the  length  of  the  upper  crack 

(region  I)  and  to  negative  perturbations  of  the  length  of  the  lower  crack  (region  III). 

For  cases  lying  in  region  II  the  system  is  stable  to  both  positive  and  negative  perturbations.  Thus,  in  this 
region  the  delaminations  will  always  grow  self-similarly.  The  solution  of  equally  spaced  delaminations  is 
represented  in  the  diagram  by  a  point  in  the  shaded  pocket  at  h^  =  h5=  1/3/z. 

In  the  case  of  unequally  spaced  delaminations,  the  constrained- contact  model  yields  a  map  with  a  larger 
region  of  stability  than  that  shown  in  Fig.  3.  The  region  is  shown  in  Fig.  D.lb.  This  difference  highlights  the 
limitations  of  the  simplified  contact  models  when  applied  to  general  delamination  geometries. 

In  a  real  structure,  the  location  of  the  delaminations  through  the  thickness  of  the  plate  will  depend  on  the 
internal  structure  of  the  material  and  the  loading  conditions,  which  therefore  determine  the  conditions  for 
the  stability  of  the  equality  of  length  of  the  system. 


Fig.  3.  Map  of  regions  of  different  energy  release  rate  for  a  system  with  two  cracks. 
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Fig.  4.  Dimensionless  diagram  of  the  critical  load  for  crack  propagation  versus  load-point  displacement  for  a  cantilever  beam  of  n 
equal  length  and  equally  spaced  cracks. 

3.3.  Macro  structural  response 

The  macro  structural  behavior  of  a  plate  of  length  10  h  with  n  equally  spaced  delaminations,  hk  = 
h/[n  +  1),  of  equal  initial  length  a0  =  2 h  is  shown  in  Fig.  4.  The  delaminations  have  been  assumed  to  prop¬ 
agate  when  fdi  =  ^cr,  where  ^cr  is  the  fracture  energy  of  the  material.  The  figure  shows  the  normalized  crit¬ 
ical  load  for  the  propagation  of  the  crack  system  as  a  function  of  the  normalized  load  point  displacement. 
As  the  number  of  delaminations  increases,  the  critical  load  for  initial  growth  of  the  delaminations  decreases 
and  the  post-peak  behavior  of  the  structure  becomes  more  ductile.  This  is  to  be  expected  as  diffuse  damage 
leads  to  a  more  ductile  structural  response  and  increased  energy  absorption.  In  contrast,  localized  damage, 
for  instance  the  propagation  of  a  single  crack,  results  in  an  increased  ultimate  capacity  of  the  structure,  but 
in  a  less  ductile  response  and  decreased  energy  absorption. 

4.  Cantilever  beam  with  two  unequal  length  delaminations 

In  a  general  system  of  delaminations  with  unequal  length  and  spacing,  Eq.  (12)  and  the  conclusions  and 
results  of  the  previous  section  do  not  apply  and  the  general  problem  defined  by  Eq.  (7)  must  be  solved  in 
order  to  define  energy  release  rates.  In  this  section,  a  system  of  only  two  cracks  is  studied,  in  which  many  of 
the  characteristics  expected  of  general  cases  can  be  observed,  without  distracting  complexity.  In  most 
cases,  the  behavior  of  the  system  depends  on  only  three  dependent  variables,  allowing  complete  visualiza¬ 
tion  of  the  interaction  effects. 

A  cantilever  beam  of  length  L  and  height  h  with  two  edge  delaminations  of  arbitrary  length  and  spacing 
is  shown  in  Fig.  5.  The  length  of  the  upper  delamination  is  au  and  the  length  of  the  lower  delamination  is 
<2L.  The  through-thickness  positions  of  the  delaminations  are  defined  by  h 3  and  h5 ,  where  h3  is  the  distance 
of  the  upper  delaminations  from  the  upper  surface  of  the  beam  and  h5  is  the  distance  of  the  lower  delam¬ 
ination  from  the  lower  surface  of  the  beam.  Utilizing  the  three  approximations  for  treating  contact,  the  sys¬ 
tem  of  equations  formed  by  the  boundary  conditions  Eqs.  (10)  and  continuity  conditions  Eqs.  (4)  together 
with  Eq.  (7)  has  been  solved  for  the  three  possible  configurations  of  the  system,  au  =  aL  (Section  3),  au  >  aL 
(Fig.  5a)  and  au  <  a L  (Fig.  5b). 

4.1.  Energy  release  rate 

Solutions  for  the  energy  release  rate  of  the  different  configurations  of  the  system  are  presented  in  Appen¬ 
dix  D.  Appendix  A  also  presents  in  Fig.  D.l  maps  constructed  for  a  general  through-thickness  distribution 
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Fig.  5.  Cantilever  beam  with  two  arbitrarily  spaced  cracks. 


and  lengths  of  the  delaminations  that  show  regions  in  which  the  energy  release  rate  of  one  of  the  cracks  is 
higher  or  lower  than  that  of  the  other  crack  and  can  be  used  to  analyze  crack  propagation  in  the  system. 

Fig.  6  shows  exemplary  solutions  for  a  system  of  two  equally  spaced  cracks  (/z3  =  h5  =  h/ 3).  Fig.  6a  de¬ 
fines  the  energy  release  rate  of  the  upper  crack  on  varying  its  length,  while  the  length  of  the  lower  crack  is 
kept  fixed  at  a^/h  =  5.0.  Fig.  6b  defines  the  energy  release  rate  of  the  lower  crack  on  varying  its  length, 
while  the  length  of  the  upper  crack  is  kept  fixed  at  a^/h  =  5.0.  The  figures  show  a  comparison  of  the  energy 
release  rates  calculated  using  the  three  approximations  of  contact  (solid  curve  =  spring-contact ,  dashed 
curve  =  unconstrained-contact  and  dash-dot  curve  =  constrained-contact ).  Numerical  results  from  these  fig¬ 
ures  and  others  show  that  the  unconstrained-and  constrained- contact  models  define  upper  and  lower  bounds 
of  the  spring-contact  model  for  all  through-thickness  distributions  of  the  cracks.  Fig.  6  also  shows  by  the 


Fig.  6.  Normalized  energy  release  rates  of  two  equally  spaced  cracks  in  the  beam  of  Fig.  5  as  a  function  of  the  length  of  the  cracks, 
(a)  Upper  crack,  (b)  lower  crack. 
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thin  line  the  energy  release  rate  of  the  crack  in  the  absence  of  the  other  crack,  or  the  single  crack  limit.  The 
limit  is  approached  when  the  interaction  effects  disappear.  The  constrained-contact  model  reaches  this  limit 
immediately  after  the  discontinuity  in  energy  release  rate.  Zheng  and  Sun  (1998)  showed  similar  behavior 
for  the  longer  delamination  of  two  central,  through-width  delaminations  in  a  three-point  bending  specimen. 

Fig.  6a  and  b  show  that  in  the  equally  spaced  delamination  system  studied  there  is  a  strong  interaction 
effect  on  the  energy  release  rate,  which  is  always  amplified  with  respect  to  the  single  crack  solution.  The  next 
section  will  show  that  different  solutions,  characterized  by  shielding,  amplification  or  a  combination  of 
both,  are  found  for  other  crack  systems.  In  addition  seems  out  of  place  reach  the  same  length  there  are  neg¬ 
ative  discontinuities  of  the  same  magnitude  in  their  energy  release  rates.  It  will  be  shown  in  Section  5  that 
the  discontinuity  is  a  good  approximation  of  a  sharp  change  in  energy  release  rate  shown  by  a  full  numer¬ 
ical  solution  of  the  problem.  In  other  material  systems  a  similar  behavior  has  been  observed.  For  instance, 
from  the  2-D  solution  of  an  infinite  body  with  a  main  crack  interacting  with  ordered  arrays  of  micro-cracks, 
sharp  jumps  in  the  stress  intensity  factors  have  been  noted  as  the  main  crack  tip  moves  through  the  array 
(Brencich  and  Carpinteri,  1996). 

4.2.  Shielding  and  amplification  of  the  energy  release  rate 

As  observed  in  the  previous  section,  a  delamination  in  a  system  of  delaminations  can  either  amplify  or 
shield  the  energy  release  rates  of  the  other  delaminations  or  have  no  influence  on  them.  This  effect  can  be 
synthesized  for  the  crack  system  of  Fig.  5  by  considering  diagrams  of  the  energy  release  rate  of  a  crack, 
normalized  with  respect  to  the  energy  release  rate  of  the  same  crack  when  the  other  crack  is  not  present,  f?io. 

The  energy  release  rate  of  a  crack  is  always  amplified  by  the  presence  of  a  shorter  crack  independent  of 
their  through-thickness  positions,  which  only  affect  the  magnitude  of  the  amplification  (the  constrained-con- 
tact  model  erroneously  predicts  for  this  regime  neither  amplification  or  shielding).  A  longer  crack  can  have 
different  effects  on  a  shorter  crack,  with  the  sense  of  the  effect  depending  on  their  through-thickness  posi¬ 
tions  and,  in  some  regimes,  their  lengths.  The  maps  shown  in  Fig.  7a  and  b  have  been  obtained  using  the 
spring-contact  model  (identical  results  are  obtained  assuming  unconstrained-contact).  The  dashed  lines  in 
the  figures  are  solutions  of  the  constrained- contact  model. 

If  the  point  describing  the  positions  of  the  two  delaminations  falls  in  the  upper  region  (III),  the  shorter 
delamination  will  always  be  shielded  by  the  longer.  The  magnitude  of  the  shielding  effect  depends  on  their 


Fig.  7.  Map  of  shielding  and  amplification  of  the  energy  release  rate  for  a  crack  in  the  presence  of  a  longer  crack  (schematic  of  Fig.  5). 
Shaded  regions  (II)  indicate  a  combination  of  shielding  and  amplification  depending  on  the  length  of  the  delaminations,  (a)  Upper 
crack,  (b)  lower  crack. 
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lengths.  If  the  point  falls  in  the  lower  region  (I),  the  energy  release  rate  of  the  shorter  delamination  will  al¬ 
ways  be  amplified  by  the  longer.  If  the  point  falls  in  the  middle  region  (II),  either  amplification  or  shielding 
can  occur  for  the  shorter  crack,  depending  on  the  two  lengths.  The  solution  of  the  constrained- contact 
model  does  not  have  a  region  of  mixed  shielding  and  amplification.  The  dashed  line  in  Fig.  7  corresponds 
to  crack  configurations  where  =  ^lo  and  separates  regions  of  shielding  (III)  and  amplification  (I).  This 
simplified  model  is  unable  to  predict  the  complex  details  of  multiply  delaminated  systems. 

An  example  considering  the  energy  release  rate  of  the  lower  delamination  for  a  through-thickness  posi¬ 
tion  of  the  cracks  that  falls  in  region  (I),  showing  amplification  of  ^L,  is  presented  in  Fig.  8a.  The  figure 
shows  normalized  with  respect  to  ^Lo  on  varying  with  the  length  of  the  upper  crack  fixed  at 
au/h  =  5.0.  When  the  lower  crack  is  shorter  than  the  upper  crack,  there  is  a  strong  amplification  of  with 


Fig.  8.  Examples  of  shielding  and  amplification  of  the  energy  release  rate  of  the  lower  crack  in  the  two-crack  system  of  Fig.  5  (upper 
crack  length  fixed). 
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an  amplification  factor  around  3.  The  amplification  is  then  reduced  by  the  discontinuity  (localized  shielding 
effect),  when  the  two  cracks  have  the  same  length,  and  it  decreases  on  increasing  the  length  of  the  crack.  An 
example  of  shielding  of  ^L?  a  point  in  region  (III),  is  shown  in  Fig.  8c.  In  this  case  the  discontinuity  in  ^ 
corresponds  to  a  localized  amplification  effect.  An  example  of  both  shielding  and  amplification  depending 
on  the  lengths  of  the  delaminations,  a  point  in  region  (II),  is  shown  in  Fig.  8b. 

In  this  structure  the  shielding  of  one  of  the  delaminations  is  always  accompanied  by  the  amplification  of 
the  other  delamination.  Among  all  the  possible  through-thickness  distributions  and  lengths  of  the  delam¬ 
inations,  the  effect  of  the  interaction  can  be  very  strong,  and  it  has  been  verified  numerically  that  the  energy 
release  rate  can  be  amplified  up  to  600%  by  the  presence  of  a  much  longer  delamination. 

Brencich  and  Carpinteri  (1996)  showed  qualitatively  similar  amplification  and  shielding  effects  when 
considering  the  interactions  of  a  main  crack  propagating  through  a  pair  of  symmetrically  located  micro¬ 
cracks  in  an  infinite  body.  For  this  problem  the  magnitude  of  the  shielding  and  amplification  depends 
on  the  spacing  between  the  cracks,  but  the  qualitative  behavior  of  the  shielding  and  amplification  is  unaf¬ 
fected  by  this  distance.  In  contrast,  for  the  structural  delaminations  discussed  in  this  paper,  due  to  geomet¬ 
rical  effects,  the  through-thickness  position  of  the  delamination  may  alter  the  condition  of  shielding  or 
amplification,  not  just  their  magnitudes.  This  is  in  agreement  with  the  numerical  study  of  the  interaction 
of  two  central  delaminations  in  a  three-point  bending  specimen  of  Zheng  and  Sun  (1998). 

4.3.  Mode  ratio 

Mode  decomposition  has  been  performed  for  the  exemplary  system  of  two  equally  spaced  delaminations. 
Fig.  9  shows  the  relative  amount  of  mode  II  to  mode  I  defined  by  the  phase  angle  ¥  =  tan_1(ifn/^i)-  Fig. 
9a  defines  Wu  for  the  upper  delamination  on  varying  its  length,  while  the  length  of  the  lower  delamination 


Fig.  9.  Relative  amount  of  mode  II  to  mode  I  stress  intensity  factors  in  the  two-crack  system  of  Fig.  5. 
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is  kept  fixed  at  a^/h  =  5.0.  Fig.  9b  defines  for  the  lower  delamination  on  varying  its  length,  while  the 
length  of  the  upper  delamination  is  kept  fixed  at  a\j/h  =  5.0. 

The  figures  compare  solutions  obtained  using  the  spring- contact  model  (thick  solid  line)  with  the  uncon¬ 
strained-contact  model  (dashed  line)  and  the  solution  of  a  single  delamination  ¥0,  in  the  absence  of  the 
other  delamination  (thin  line).  For  a  single  delamination,  the  relative  amount  of  mode  II  to  mode  I  is  inde¬ 
pendent  of  the  delamination  length  and  given  by  ^uq  =  66.7°  for  the  upper  delamination  and  by  *FL0  =  90° 
for  the  lower  delamination  (the  lower  delamination  experiences  pure  mode  II  conditions  due  to  localized 
contact  at  the  crack  tip). 

Fig.  9a  shows  that  when  the  upper  delamination  is  shorter  than  the  lower,  the  conditions  are  pure  mode 
II.  Comparison  with  the  single  delamination  solution  highlights  the  strong  effect  produced  by  the  presence 
of  the  lower  delamination.  When  the  upper  delamination  reaches  the  length  of  the  lower  delamination  there 
is  a  sharp  transition  and  the  mode  I  component  is  suddenly  increased  above  the  corresponding  single 
delamination  solution  (thin  line).  As  the  upper  crack  lengthens,  the  interaction  decreases  and  the  solution 
tends  to  the  single  delamination  limit.  The  dashed  curve  highlights  the  limitations  of  the  simplified  uncon¬ 
strained-contact  model,  which  predicts  an  unrealistic  phase  angle  ¥  >  90°,  due  to  the  interpenetration  of  the 
delamination  surfaces. 

Fig.  9b  shows  that  the  lower  delamination  is  in  pure  mode  II  conditions  due  to  contact  at  the  delami¬ 
nation  tip  when  its  length  is  shorter  than  the  length  of  the  upper  delamination  and  this  solution  coincides 
with  the  solution  of  the  single  delamination.  When  the  lower  delamination  approaches  the  length  of  the 
upper  delamination,  there  is  a  sudden  transition  (solid  curve)  to  a  large  value  of  the  mode  I  component. 
The  interaction  effect  then  disappears  quickly  as  the  lower  delamination  lengthens  and  the  solution  tends 
again  to  the  solution  of  a  single  delamination.  Interestingly,  when  the  lower  delamination  is  longer,  the 
bending  theory  model  always  predicts  opening  at  the  lower  delamination  tip  even  when  the  mode  decom¬ 
position  method  defines  a  negative  mode  I  stress  intensity  factor  and  therefore  pure  mode  II  conditions 
(shown  in  Fig.  9b  for  a^/h  >6.75).  This  fictitious  opening  displacement  field  is  a  consequence  of  neglecting 
the  root  rotations  at  the  delamination  tip,  which  if  properly  accounted  for  would  restore  contact  (Andrews 
et  al.,  2005). 

4.4.  Delamination  growth  and  macro- structural  behavior 

To  investigate  the  macro-structural  response  of  the  plate,  the  quasi-static  propagation  of  the  system  of 
two  delaminations  shown  in  Fig.  5  has  been  studied  using  the  spring-contact  model.  A  delamination  has 
been  assumed  to  propagate  when  its  energy  release  rate  equals  the  fracture  energy  of  the  material,  ^cr. 
The  cracks  have  been  assumed  to  grow  under  delamination  length  control,  which  allows  virtual  delamina¬ 
tion  growth  branches  to  be  followed.  Several  cases,  identified  by  different  through-thickness  positions  of  the 
delaminations  and  different  notch  lengths  as  shown  in  Figs.  10-13,  have  been  considered  to  highlight  a 
number  of  interesting  macro  structural  responses  of  the  system.  The  figures  show  the  critical  load  for  crack 
propagation  as  a  function  of  the  load  point  deflection.  Also  shown  in  the  inset  in  the  figures  is  the  evolution 
of  the  lengths  of  the  delaminations;  the  dashed  line  shows  the  length  of  the  lower  delamination  and  the  solid 
line  shows  the  length  of  the  upper  delamination  during  the  loading  process. 

For  the  first  delamination  configuration,  shown  in  Fig.  10,  the  lower  delamination  starts  to  propagate 
first  at  (A).  The  propagation  is  unstable  up  to  point  (B),  where  the  delamination  reaches  the  upper  delam¬ 
ination.  At  this  point  the  delamination  arrests  and  the  delamination  system  can  be  made  to  propagate  only 
by  increasing  the  applied  load.  This  behavior  is  due  to  a  negative  discontinuity,  shielding,  in  the  energy  re¬ 
lease  rate.  After  point  (C)  the  lower  delamination  continues  to  propagate  unstably  until  the  structure  fails. 
A  snap-back  instability  is  present  as  the  lower  delamination  grows  to  the  upper  delamination. 

In  the  case  shown  in  Fig.  1 1  the  lower  delamination  begins  to  propagate  unstably  at  point  (A).  At  (B)  the 
delaminations  are  of  the  same  length  and  there  is  a  sudden  drop  in  the  critical  load,  corresponding  to  a 


M.G.  Andrews  et  al.  I  International  Journal  of  Solids  and  Structures  43  (2006)  853-886 


871 


A,  =  0.25A  A,  =  0.4A 
aL  =  2.5A  L  =  10h 


JoJiE 


Fig.  10.  Dimensionless  diagram  of  the  critical  load  for  crack  propagation  versus  load-point  displacement  in  the  two-crack  system  of 
Fig.  5:  snap-back  instability  and  local  increase  in  critical  load  due  to  a  local  shielding  effect. 
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Fig.  1 1 .  Dimensionless  diagram  of  the  critical  load  for  crack  propagation  versus  load-point  displacement  in  the  two-crack  system  of 
Fig.  5:  snap-back  instability  with  local  drop  in  critical  load  due  to  a  local  amplification  effect. 
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Fig.  12.  Dimensionless  diagram  of  the  critical  load  for  crack  propagation  versus  load-point  displacement  in  the  two-crack  system  of 
Fig.  5:  snap-through  instability  and  hyper-strength  phenomenon  due  to  a  local  shielding  effect. 
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Fig.  13.  Dimensionless  diagram  of  the  critical  load  for  crack  propagation  versus  load-point  displacement  in  the  two-crack  system  of 
Fig.  5:  crack  pull  along. 


positive  discontinuity,  amplification,  of  the  energy  release  rate.  The  lower  delamination  then  continues  to 
propagate  unstably. 

The  case  shown  in  Fig.  12  is  similar  to  that  examined  in  Fig.  10.  However  in  this  case  the  new  critical 
load  due  to  the  negative  discontinuity,  shielding,  in  the  energy  release  rate  is  higher  than  the  initial  load  for 
delamination  propagation.  This  is  an  interesting  example  of  hyperstrength,  which  is  defined  as  the  condi¬ 
tion  that  the  ultimate  load  of  the  structure  exceeds  the  load  corresponding  to  the  first  propagation  of  one  of 
the  delaminations  of  the  system  even  in  the  absence  of  any  strengthening  mechanisms.  The  hyperstrength  is 
produced  by  the  elastic  interactions  of  the  parallel  delaminations  that  create  local  shielding  and  hardening 
effects,  leading  to  a  snap- through  instability. 

An  example  of  crack  pull-along,  in  which  the  finite  propagation  of  one  delamination  causes  the  propa¬ 
gation  of  the  other  delamination,  is  shown  in  Fig.  13.  In  this  case  the  lower  delamination  propagates  first  at 
point  (A).  When  it  reaches  the  second  delamination  at  point  (B),  there  is  an  increase  in  critical  load  due  to  a 
negative  discontinuity  in  energy  release  rate.  The  lower  delamination  continues  to  propagate  at  point  (C). 
At  point  (D),  the  upper  delamination  begins  to  propagate,  and  both  delaminations  continue  to  propagate, 
keeping  a  constant  relative  length  until  the  structure  fails.  This  behavior  occurs  for  through-thickness  dis¬ 
tributions  of  the  delaminations  that  fall  just  to  the  right  or  just  above  the  shaded  region  of  the  map  of  Fig. 
D.la,  and  just  to  the  right  of  the  shaded  region  in  the  map  of  Fig.  D.lb. 

5.  Validation  of  the  proposed  model 

The  proposed  model,  applied  to  the  problem  of  a  cantilever  beam  with  two  edge  cracks,  has  been  val¬ 
idated  using  the  finite  element  method  (FEM)  and  the  commercial  code  ANSYS  (5.5).  The  mesh  consisted 
of  plane  stress  isoparametric  triangular  elements.  To  obtain  an  accurate  solution,  including  relative  crack 
displacements,  while  keeping  the  number  of  degrees  of  freedom  low,  a  coarse  mesh  was  used  for  the  body  of 
the  beam  and  very  fine  meshes  around  each  crack  tip  and  along  areas  of  contact.  The  stress  singularities  at 
the  crack  tips  were  modeled  with  rosettes  of  quarter  point  elements.  Contact  along  the  crack  faces  was  sim¬ 
ulated  using  gap  elements  that  prevent  interpenetration  of  the  beams  with  stiff  linear  springs.  Convergence 
of  the  solution  was  checked  by  varying  the  size  and  number  of  elements  and  stiffnesses  of  the  gap  springs. 
Energy  release  rates  were  calculated  by  two  methods,  the  /-Integral,  Rice  (1968),  and  the  displacement  cor¬ 
relation  method,  Chan  et  al.  (1970),  which  was  also  used  to  evaluate  the  mode  ratio.  The  path  of  the  /-inte¬ 
gral  was  chosen  within  a  fine  mesh  region  so  that  it  encompassed  sufficient  elements  for  its  value  to 
converge. 
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5.1.  Comparison  of  energy  release  rates 

The  energy  release  rates  determined  by  the  finite  element  method  and  the  present  spring-contact  model 
are  shown  in  Fig.  14a  for  a  system  of  two  equally  spaced  cracks  (/z3  =  h5  =  /z/3).  The  error  between  the  two 
solutions  is  shown  in  Fig.  14b.  In  the  figure,  the  normalized  energy  release  rate  of  the  upper  crack  is 
depicted  as  a  function  of  the  length  of  the  upper  crack.  The  length  of  the  lower  crack  is  kept  fixed  at 
a^/h  =  5.0. 

Fig.  14a  and  b  show  that  the  results  of  the  proposed  model  generally  agree  with  the  finite  element  results. 
When  the  cracks  reach  the  same  length,  the  finite  element  solution  predicts  a  continuous  transition  in  en¬ 
ergy  release  rate.  This  transition  matches  the  discontinuity  in  energy  release  rate  predicted  by  the  beam  the¬ 
ory  model. 

Short  cracks  (a\j  <  2 h  for  the  geometry  of  Fig.  14)  are  characterized  by  large  errors  due  to  the  well- 
known  limitations  of  classical  beam  theory.  As  expected,  the  error  drops  when  the  length  of  the  crack  in¬ 
creases,  provided  it  is  not  similar  to  the  length  of  the  other  crack,  to  less  than  5%  when  2h<  a\j<  4.6/z  or 
flu  >  7.0/l  Where  the  two  cracks  have  similar  lengths  and  the  energy  release  rate  is  characterized  by  a  dis¬ 
continuity,  the  error  increases  significantly  (up  to  30%  for  the  geometry  studied  in  Fig.  14).  The  width  of 
this  region  is  given  by  a  few  times  the  separation  of  the  planes  of  the  interacting  cracks  and  depends  on  the 
conditions  of  the  crack  surfaces  (contact  or  opening).  For  the  geometry  studied  in  Fig.  14,  the  error  is  above 
5%  for  4.6/z  <  a\j  <1  .Oh.  The  interval  is  not  symmetric  about  the  length  of  the  lower  crack,  aL  =  5/z,  because 
of  the  different  conditions  along  the  crack  surfaces  when  the  upper  crack  is  shorter  (opening)  or  longer 
(contact)  than  the  lower.  Numerical  calculations  show  that  behaviors  similar  to  that  observed  in  Fig.  14 
characterize  all  possible  crack  configurations  and  that  the  width  of  the  region  where  the  error  increases  sig¬ 
nificantly  is  always  a  few  times  the  distance  between  the  planes  of  the  cracks. 
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Fig.  14.  (a)  Energy  release  rate  of  the  upper  crack  in  the  two-crack  system  of  Fig.  5:  comparison  between  FE  and  proposed  model 
results,  (b)  Relative  error  between  the  solutions  of  (a)  (dashed  line  indicates  a  5%  relative  error). 
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5.2.  Comparison  of  mode  ratios 

The  relative  amount  of  mode  II  to  mode  I  as  measured  by  the  phase  angle  W  =  tan-1  (Kn/Kj)  determined 
by  the  finite  element  method  and  the  model  presented  in  this  paper  utilizing  the  spring- contact  approxima¬ 
tion  is  shown  in  Fig.  15.  Fig.  15a  shows  the  phase  angle  of  the  upper  crack  as  a  function  of  the  length  of  the 
upper  crack,  while  keeping  the  length  of  the  lower  crack  fixed  at  ajh  =  5.0.  Fig.  15b  shows  the  phase  angle 
of  the  lower  crack  while  keeping  the  length  of  the  upper  crack  fixed  at  a^/h  =  5.0. 

The  phase  angles  of  the  lower  crack  determined  by  the  two  models,  shown  in  Fig.  15b,  agree  well  with 
the  relative  error  peaking  near  the  discontinuity  at  only  10%.  Fig.  15a  instead  shows  a  significant  difference, 
with  an  error  up  to  26%,  between  the  results  of  the  bending  theory  model  and  FEM  when  the  upper  crack  is 
longer  and  there  is  contact  along  both  delamination  surfaces.  When  the  lower  crack  is  longer  and  there  is  no 
contact  along  the  delamination  surfaces,  the  error  is  instead  very  low.  This  result  confirms  what  has  already 
been  observed  for  the  energy  release  rate.  However,  the  width  of  the  region  where  the  relative  error  on  the 
mode  ratio  is  quite  large  and  is  greater  than  that  corresponding  to  the  energy  release  rate. 

5.3.  Discussion 

The  relative  error  between  the  model  and  the  FE  predictions,  as  observed  in  Figs.  14  and  15,  is  greater 
when  there  is  contact  between  the  crack  surfaces.  This  is  mainly  due  to  the  assumption  of  the  proposed 
model  of  zero  relative  root  rotations  of  the  sublaminates  at  the  crack  tips  that  leads  to  an  overestimation 
of  the  contact  pressures.  The  effects  of  the  root  rotations  on  the  contact  pressures  and  thus  on  the  energy 
release  rate  and  mode  ratio  are  most  significant  when  the  cracks  are  close  to  the  same  length  \ajj  —  <zL|  <  h 
and  are  stronger  on  the  mode  ratio.  When  there  is  no  contact  along  the  crack  faces,  neglecting  the  root 


Fig.  15.  Relative  amount  of  mode  II  to  mode  I  for  a  two-crack  system:  comparison  between  proposed  model  and  FE  results. 
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Fig.  16.  Comparison  of  predicted  contact  pressures  along  the  faces  of  the  upper  crack  in  the  two-crack  system  (dashed  line:  contact 
pressures  using  the  FE  displacement  field  and  spring-contact  approximation). 


rotations  does  not  significantly  affect  the  solution,  confirming  what  has  already  been  observed  for  single 
delaminations  loaded  primarily  in  shear. 

The  contact  pressure  distribution  along  the  surfaces  of  the  upper  crack,  for  equally  spaced  cracks  with 
lengths  ajj/h  =  5.5  and  a^/h  =  5,  predicted  by  the  proposed  model  and  the  FE  method  are  shown  in  Fig.  16 
as  a  function  of  the  position  along  the  crack.  The  model  predicts  a  peak  pressure  double  that  of  the  FE 
solution.  These  discrepancies  in  the  solutions  cannot  be  attributed  to  the  contact  approximation  used  in 
the  model:  using  the  centerline  beam  deflections  obtained  in  the  FE  solution  and  Eq.  (6)  to  predict  the  con¬ 
tact  pressure  leads  to  the  pressures  identified  by  the  dashed  curve  that  for  magnitude  and  shape  closely  fol¬ 
low  the  finite  element  solution.  Including  shear  deformations  in  the  model  improves  the  pressure 
distribution,  but  does  not  substantially  improve  the  values  of  the  energy  release  rate. 

The  overestimation  of  the  contact  pressures  is  associated  with  significant  errors  in  the  stress  resultants  at 
the  crack  tip.  This  is  shown  in  Table  1  where  bending  moment  and  shear  force  obtained  using  the  proposed 
model  are  compared  with  the  FE  results  (rows  2  and  1).  Table  1  also  compares  the  energy  release  rates  and 
phase  angles  obtained  using  the  proposed  model  (row  2)  and  the  FEM  (row  1)  and  highlights  the  previously 
noted  large  errors.  If  the  actual  FE  stress  resultants  are  used  instead  of  the  model  resultants  in  Eq.  (9a)  and 
in  the  expressions  of  Suo  and  Hutchinson  for  the  stress  intensity  factors,  the  solution  is  still  affected  by  large 
errors.  This  is  shown  in  row  3  of  the  table.  The  reason  for  this  discrepancy  is  that  both  Eq.  (9a)  for  the 
energy  release  rate  and  Suo  and  Hutchinson’s  expressions  do  not  account  for  the  contributions  due  to 
the  shear  deformations  along  the  beams  and  the  deformations  at  the  crack  tip  cross  sections  (i.e.,  root 
rotations). 

The  expression  (9a)  for  the  energy  release  rate  can  be  modified  to  account  for  these  effects  as  follows: 


M) 

Eli 


V2  N2  N 

— - - J-  +  2Vicpi 

kGAj  EAj  J^J 


(13) 


where  ^  is  the  shear  modulus,  k  is  the  shear  correction  coefficient  (5/6  for  a  rectangular  cross  section), 
V  and  (p  are  the  shear  resultant  and  the  bending  rotation  respectively. 

Li  et  al.  (2004)  recently  modified  the  stress  intensity  factor  expressions  derived  by  Suo  and  Hutchinson  to 
include  the  contributions  of  the  shear  deformations  along  the  beam  and  the  crack  tip  deformations  numer¬ 
ically.  Wang  and  Qiao  (2004)  proposed  an  analytical  extension  of  the  method  of  Suo  and  Hutchinson. 
Rows  (4)  and  (5)  in  the  Table  show  solutions  obtained  using  these  two  methods.  Both  methods  lead  to  sub¬ 
stantial  improvement  in  the  solutions.  The  method  of  Wand  and  Qiao  relies  on  the  assumptions  of  plate 
theory  in  order  to  determine  the  crack  tip  rotations,  and  thus  is  not  as  accurate  as  the  numerical  solution 
of  Li  et  al. 
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Table  1 

Stress  resultants,  energy  release  rate  and  phase  angle  of  the  upper  crack 


Moment — 

Moment — 

Shear — 

Shear — 

^u, 

gvEh/P2, 

upper 

lower 

upper 

lower 

%  Error  from 

%  Error  from  FEM 

beam  arm 

beam  arm 

beam  arm 

beam  arm 

FEM 

(M/Ph) 

(M/Ph) 

(V/P) 

(V/P) 

(1) 

FE  stress  resultants 
^  and  W  Displacement 
correlation  Technique 

-0.967 

-4.535 

-0.644 

1.644 

72.4° 

428.59 

(2) 

Proposed  model 

-0.300 

-5.200 

-1.752 

2.752 

54.2°, 

380.89, 

25.18% 

11.13% 

(3) 

FE  stress  resultants 

-0.967 

-4.535 

-0.644 

1.644 

80.9°, 

386.14, 

^  Eq.  (9a) 

W  after  Suo  and  Hutchinson 

-11.66% 

9.91% 

(4) 

FE  stress  resultants 

-0.967 

-4.535 

-0.644 

1.644 

72.2°, 

430.18, 

W  after  Li  and  Thouless 

0.31% 

-0.37% 

(5) 

FE  stress  resultants 

-0.967 

-4.535 

-0.644 

1.644 

74.4°, 

417.65, 

0  Eq.  (13) 

W  after  Wang  and  Qiao 

-2.71% 

2.55% 

Equally  spaced  cracks,  au/h  =  5.5,  ajh  =  5,  h/L  =  0.1. 


The  observations  above  highlight  that  the  model  proposed  in  this  paper  could  be  improved  substantially 
by:  removing  the  assumption  of  zero  relative  root  rotations  at  the  crack  tips,  for  instance  by  using  localized 
linear  elastic  rotational  springs  (Pandey  and  Sun,  1996);  including  the  contribution  of  shear  deformations 
along  the  beams;  applying  Eq.  (13)  for  the  energy  release  rate  and  an  extended  solution,  e.g.  of  Li  et  al. 
(2004)  or  Wang  and  Qiao  (2004),  for  the  decomposition  of  the  modes  of  fracture.  This  would  lead  to  a  more 
accurate  quantitative  prediction  of  the  displacement  fields  along  the  crack  surfaces  and  the  fracture  param¬ 
eters  of  the  system.  Such  an  approach  would  become  necessary,  and  would  probably  be  sufficient  for  plane 
conditions,  if  the  influence  of  crack-wake  mechanisms,  such  as  bridging  or  cohesive  mechanisms  and  fric¬ 
tion,  on  the  fracture  behavior  of  multiply  delaminated  plates  is  to  be  investigated. 


6.  Application  of  the  proposed  model  to  the  design 

6.1.  Beam  theory  for  quantitative  analyses 

Elastic  analyses  of  laminates  that  are  based  on  plate  elements  are  widely  used  and  accurate  for  many 
design  predictions.  Extended  laminar  sheets  comprising  thin  plies  do  indeed  very  nearly  satisfy  the  condi¬ 
tions  of  strain  distribution  necessary  for  the  reduced  degrees  of  freedom  in  beam  (or  plate)  theory  to  be 
sufficient. 

As  shown  in  this  paper,  modeling  delamination  fracture  is,  however,  much  more  challenging.  In  multiply 
delaminated  beams  as  well  as  other  more  complex  structures,  delaminations  propagate  in  mixed  mode,  of¬ 
ten  predominantly  in  shear  (modes  II  and  III).  Crack  propagation,  even  in  the  absence  of  a  through-thick¬ 
ness  reinforcement  (stitching,  pins,  etc.),  will  therefore  be  very  strongly  and  commonly  influenced  by 
contact  forces  and  crack  face  friction.  Other  crack  wake  effects,  such  as  bridging  by  fibers  or  craze  fibrils, 
might  also  be  important.  The  tractions  imposed  on  the  fracture  surfaces  by  such  effects  will  depend,  pos¬ 
sibly  sensitively,  on  the  mixed  mode  crack  displacement  vector;  and  therefore  accurate  prediction  of  the 
crack  displacement  becomes  a  prerequisite  to  a  full  quantitative  model  of  delamination. 
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The  question  of  whether  plate  or  beam  theory  methods  can  deal  with  quantitative  delamination  analysis 
turns  on  the  success  of  estimates  of  root  rotations,  which  control  the  accuracy  of  crack  displacement  predic¬ 
tions  near  crack  tips.  Our  current  work  and  work  in  the  literature  (Andrews  et  al.,  2005;  Sun  and  Pandey, 
1994;  Pandey  and  Sun,  1996)  encourages  the  view  that  reasonably  simple  root  rotation  corrections  can  in¬ 
deed  restore  accuracy.  If  so,  beam  theory  may  be  not  only  useful  but  even  the  preferred  method,  especially 
for  plane  problems,  including  many  common  experimental  tests,  since  it  yields  straightforward,  relatively 
small  numerical  problems  and,  in  many  important  cases  and  limits,  highly  instructive  analytical  results. 

Competition  in  accurate  delamination  modeling  will  come  from  modern  cohesive  element  methods, 
which  are  extensions  of  finite  element  formulations  to  include  elements  that  can  introduce  bridging  and 
contact  tractions  in  convenient  ways.  Cohesive  models  are  an  attractive,  physically  based  formulation  of 
nonlinear  crack  tip  and  crack  wake  processes,  very  well  suited  to  complex  3D  configurations  with  arbitrary 
mode  mixtures.  They  have  a  growing  record  of  successful  simulations  of  experiments  (e.g.,  Remmers  et  al., 
2003;  Yang  and  Cox,  2004).  However,  accurate  cohesive  model  analyses  of  complex  laminated  structures 
may  comprise  106  degrees  of  freedom  and  do  not  yield  analytical  results,  even  in  limits.  Beam  theory 
and  computationally  intensive  cohesive  element  methods  may  prove  complementary. 

6.2.  Towards  structural  design  concepts 

Many  configurations  of  loads  and  boundary  conditions  arise  in  structures  other  than  the  cantilever  case 
studied  here.  Nevertheless,  the  cantilever  case  suggests  some  principles  concerning  the  role  of  the  multiplic¬ 
ity  of  delaminations  in  structural  performance  that  may  prove  generally  applicable  to  cases  of  mixed  shear 
and  bending. 

Assume  that  residual  strength  following  delamination  damage  declines  as  the  size  of  the  cracks  increases 
(e.g.  Fig.  4).  Then  the  results  for  uniformly  spaced  cracks  (Eq.  (12))  suggest  that,  to  enhance  residual 
strength,  material  conditions  that  favour  greater  number  of  cracks  uniformly  distributed  through  a  structure 
are  to  be  avoided.  On  the  other  hand,  if  the  cracks  are  confined  to  a  band  of  limited  width,  then  resistance  to 
growth  rises  with  their  number.  Protection  against  long  cracks  might  therefore  be  sought  by  deliberately 
introducing  multiple  weaker  planes  within  a  limited  band  in  the  material.  Of  course,  the  effect  of  such  a  strat¬ 
egy  on  other  failure  routes  must  also  be  considered.  For  instance  by  tailoring  the  material  for  increased  resid¬ 
ual  strength  the  ductility  of  the  structure  is  reduced,  leading  to  a  more  brittle  failure  of  the  structure. 

The  case  study  of  the  system  of  two  cracks  shows  that  interaction  effects  are  strong.  Cases  in  which  one 
crack  shields  the  other  and  cases  in  which  it  leads  to  accelerated  growth  in  its  partner  can  be  found,  depend¬ 
ing  on  where  the  two  cracks  reside  in  the  laminate.  The  possibility  of  acceleration  implies  that  design  and 
life  predictions  based  on  solutions  for  a  single  crack  cannot  be  safe. 

In  many  cases  studied  here,  the  possibility  of  crack  face  contact  arises.  Long  zones  of  crack  wake  contact 
and  friction  can  therefore  be  expected  to  be  common  in  structural  problems  involving  shear  and  bending. 
The  friction  zones  will  not  generally  be  confined  to  a  near-tip  region,  in  particular  one  that  is  sufficiently 
small  that  friction  effects  could  be  incorporated  by  modifying  the  critical  value  of  the  energy  release  rate. 
Rather,  the  friction  contact  zone  will  vary  in  extent  (and  possibly  in  the  number  of  contact  domains)  with 
the  crack  size  and  its  relation  to  other  cracks.  The  contact  zones  must  be  calculated  explicitly.  Similar  re¬ 
marks  will  apply  to  treating  through- thickness  reinforcement,  should  it  be  present.  Work  is  in  progress  to 
investigate  these  topics. 


7.  Conclusions 

A  model  has  been  presented  that  allows  for  the  analysis  of  laminated  plates  with  multiple  delaminations 
deforming  in  cylindrical  bending.  The  model  utilizes  the  assumptions  of  classical  beam  theory,  which 
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neglects  shear  and  near-tip  deformation,  and  accounts  for  non-frictional  contact  between  the  crack  sur¬ 
faces.  It  has  been  applied  to  a  cantilever  beam  with  multiple  through-width  edge  delaminations  subject 
to  a  concentrated  load  at  its  free  end.  Contact  has  been  simulated  using  different  approaches,  uncon¬ 
strained -,  constrained-  and  spring-contact ,  where  the  names  reflect  the  approximation  used  to  describe  the 
displacement  field.  Finite  element  analyses  applied  to  a  two  delamination  system  show  that  the  spring-con- 
tact  model  is  the  most  accurate;  however  utilizing  it  results  in  a  non-linear  problem  and  thus  requires  a  par¬ 
tially  numeric  solution  while  the  other  models  result  in  closed  form  solutions.  The  two  simpler 
approximations  give  upper  and  lower  bounds  of  the  energy  release  rate  determined  by  the  spring- contact 
approximation.  In  the  case  of  the  mode  ratio,  however  they  give  results  that  are  often  qualitatively  different 
and  at  times  also  incorrect.  These  problems  are  expected  to  be  compounded  in  systems  of  more  delamin¬ 
ations  with  more  complicated  states  of  contact.  Thus  when  accurate  results  are  required,  the  simplified  con¬ 
tact  approximations  should  only  be  used  as  general  bounds  and  as  a  first  approximation  of  the  regions  of 
contact  for  the  iterative  solution  of  the  spring-contact  model. 

Analysis  of  the  energy  release  rate  and  mode  ratio  of  the  delaminations  has  revealed  several  key  insights 
into  the  behavior  of  multiply  delaminated  systems.  When  delamination  growth  is  considered,  a  discontinu¬ 
ity  in  the  energy  release  rate  and  mode  ratio  is  found  when  the  delaminations  reach  the  same  length.  This 
discontinuity  leads  to  instantaneous  shielding  or  amplification  of  the  energy  release  rate,  eliciting  behaviors 
such  as  local  snap-back  and  snap-through  instabilities,  hyperstrength,  crack  pull-along  and  crack  arrest.  All 
of  these  behaviors  are  strongly  dependent  on  and  controlled  by  the  through-thickness  spacing  of  the 
delaminations.  The  energy  release  rates  of  the  delaminations  are  also  amplified  or  shielded  when  the  delam¬ 
inations  are  not  of  the  same  length,  again  depending  on  the  through-thickness  positions  of  the 
delaminations.  These  effects  show  some  similarity  to  the  effects  previously  observed  in  the  interaction  of 
microcracks  or  of  main  cracks  with  clouds  of  microcracks  in  infinite  media,  which  can  represent  damage 
in  concrete  and  coarse  grain  ceramics  (Hutchinson,  1987;  Rose,  1986;  Rubinstein,  1985;  Kachanov, 
1986;  Brencich  and  Carpinteri,  1996).  However  in  the  problem  considered  here,  strong  modifications  in 
the  results  and  more  complex  behaviors  are  observed  due  to  the  finiteness  of  the  structure. 

The  through-thickness  positions  of  the  delaminations  also  control  the  sensitivity  of  a  system  of  equal 
length  delaminations  to  length  perturbations.  When  the  delaminations  are  equally  spaced  the  equality  of 
length  is  stable  (i.e.  simultaneous  growth  of  the  delaminations)  leading  to  a  more  ductile  failure;  when  they 
are  not,  the  equality  of  length  is  often,  but  not  always,  unstable  resulting  in  a  more  brittle  response  (i.e. 
growth  of  one  or  a  limited  number  of  delaminations). 

While  the  present  study  was  limited  to  a  system  of  only  two  unequal  length  delaminations,  the  proposed 
model  could  be  applied  to  analyze  systems  with  a  general  number  of  delaminations  following  the  procedure 
described  in  Section  2.  Behaviors  similar  to  those  of  the  two-delamination  system  are  expected  in  a  system 
with  many  delaminations,  resulting  in  a  macro  structural  response  with  a  saw-tooth  appearance,  with  re¬ 
peated  increases  and  decreases  in  critical  load  as  the  critical  delamination  tip  grows  with  and  past  other 
delamination  tips  in  the  system.  Different  boundary  and  loading  conditions  from  those  studied  here  are  ex¬ 
pected  to  have  a  mostly  quantitative  effect  on  the  results  and  can  be  easily  solved  using  the  same  model  and 
solution  procedure.  With  minimal  modification,  friction  acting  in  the  regions  of  contact  can  also  be  con¬ 
sidered.  Its  effect  will  be  additional  shielding  of  the  energy  release  rate  at  the  delamination  tips. 
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Appendix  A.  Derivation  of  general  solutions  for  system  of  coupled  differential  equations 

In  this  appendix,  the  general  solutions  for  the  system  of  differential  equations,  Eqs.  (7)  and  (8)  are  de¬ 
rived,  for  one  free  segment  and  for  2  and  3  beam  segments  in  contact. 

A.l.  Single  beam  segment 


For  a  single  beam  segment  A  the  governing  differential  equations,  Eqs.  (7)  and  (8),  are 

EIaWa^xxx  =  0, 

EAaua,xx  =  0 

the  general  solutions  of  this  system  of  equation  is 

Wa  (x)  =  C\  C2X  ~b  ~\~  C^x^, 
uA  (x)  =  C5  +  C6x, 

where  the  Cs  are  unknown  constants  of  integration. 


(A.l) 


(A.2) 


A. 2.  Two  beam  segments  in  contact 


For  two  beam  segments  A  and  B  in  contact,  the  governing  differential  equations,  Eqs.  (7)  and  (8)  are 


EIa^Axxxx  +  kA}B(wA  —  wB)  —  0, 
EAaua,xx  —  0, 

EIBwBxxxx  —  kA)B(wA  —  wB)  =  0, 
EABuB:XX  =  0, 


(A.3) 


where  kA  B  is  given  by  Eq.  (5).  Solutions  of  the  form: 

wA(x)  =A,eKX,  wB(x)  =  A2  eKX  (A.4) 

are  sought.  The  resulting  algebraic  equation,  determined  by  substitution  of  Eq.  (A.4)  into  Eqs.  (A.3)  and 
subsequent  elimination  of  A 1  or  A2  is 

*V  +  4,Qe“  =  0,  (A.5) 

where 


The  relationship  between  the  unknown  constants  A\  and  A2  is 

A2  EIb  4 

T  =  1 — K  +l- 
A 1  kA,B 

The  8  roots  of  Eq.  (A.5)  are 

K  =  0,0,0,0,±(1 


(A.6) 


(A.7) 

(A.8) 
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where  i  =  yf—Y.  The  general  solution  of  this  system  of  differential  equations  is  therefore 

wA(x)  =  El A  e^x(Ci  cos  PABx  +  C2  sin  fiABx) 

H-  El A  e  ^a’bX(C3  cos  fiABx  H-  C4  sin  fiA  Bx^)  H-  C$x^  H-  C^x^  -f  C7X  -f  C8,  (A. 9) 

uA{x)  =  C9  +  Ciox 

w5(x)  =  e^’fiX(Ci  cos  fiABx  +  C2  sin  PAtBx), 

H-  EIbc  ^ ‘a,bX^(J 3  qos  fiABx  H-  C4  sin  fiA  Bx^  H-  C$x^  -f  +  C7X  +  C8  (A. 10) 

uA{x)  =  C 11  +  Ci2x. 


A3.  TTzree  segments  in  contact 


The  governing  differential  equations  for  three  beam  segments  v4,  B  and  C  in  contact  are 

El  A  wAjXXXX  +  kAB  (wA  —  wB)  =  0, 

EAauAxx  0, 

EIbWB,xxxx  +  kB,c{wB  —  wc)  =  kA)B(wA  —  wB ), 

EABUb.xx  =  0, 

EICWCjcxxx  —  kB,c{wB  —  wc)  =  0, 


EAcuC:xx  =  0. 

The  /:’s  are  given  by  Eq.  (A.5).  Solutions  of  the  form 
wA(x)  =  A\ eKX,  w5(x)  =  A2eKX,  wc(x)  =A3eKX 
are  sought.  The  resulting  algebraic  equation  for  determination  of  k  is 

k4(k8  +  a a,b,ck4  +  Vaac)  e“  =  °> 

where 

OCAAC  =  kA,B  (-ETT  +  +  kB,C  +  AjE\ , 


yA,B,c  —  kA,BkB:c 


El A  EIb 

1 


-  +  - 


1 


-  +  - 


1 


\EIaEIb  EIaEIc  E1bEICj 
The  relationships  between  the  unknown  constants  A\,  A2  and  A 3  are 


^2  EIC  4 
— -  =  7 - Tv  T  1 , 


A 


^B,C 


A 1  EIbEIc  8 

—  = - k; 

A  kABkB  c 


-'-(EIb+EIc) 

Ka,C 


-  El  r  TC 


The  12  roots  of  Eq.  (A.  13)  are 

3 


k  =  0,0,0,0,  ±(1  ±  i)  (  -uA,B,c  ±- 


'y4,5,C 


-4y 


^,5,C 


1/4 


(A.ll) 


(A.12) 
(A.  13) 


(A.  14) 


(A.  15) 


(A- 16) 
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The  general  solutions  of  the  system  of  equations  is  therefore: 

wA  (x)  =  C\  T  C2x  C3x^  T  C/pc1 

+  p4e  ^ X(C5  sin  ffp^x  +  C6  cos  ,yp^x)  +  p3  e  ^*(C7  sin  ^yp^x  +  C8  cos  .J/p^x) 

+  p3  e“  ^ (C9  sin  ffp[x  +  Ci0  cos  ffp~xx)  +  p4 e“  ^x(Cn  sin  .yp^x  +  Ci2  cos  </p^x) , 
Wb{x)  —  C*i  - h  C2X  -|-  C3x^  T  C4X^ 

+  p6 e  ^X(C5  sin  ^yp^x  +  C6  cos  ^yp»  +  p5  e  ^X(C7  sin  ^y^x  +  C8  cos  ^J/^x) 

+  p5e“  ^X(C9  sin  ^/p^x  +  Cxo  cos  ^yp^x)  +  p6e_  sin  ^yp^x  +  Cu  cos  ,yp^x), 

Wc(x)  =  C\  T  C2x  T-  C3x^  T  C*4x^ 

+  e  ^X(C5  sin  -^/p^x  +  C6  cos  ffppx)  +  e  ^X(C7  sin  ffp~xx  +  C8  cos  ^/ppc) 

+  e“  ^X(C9  sin  ^yp^x  +  Cxo  cos  ^/p^x)  +  e-  ^x(Cn  sin  ^/p^x  +  Cn  cos  yffcx), 

where 


Pi  5  P2  —  g  a^,5,C  T  \j ^A,B,C  ^^A,B^B,Cy A^B^C'i 

16  EIBEIcp\  El BkBC  +  El  ckBC  +  El  ckA,B 

P3  =  - L-4 - 1 - - “Pi  +  1, 


P4  = 


^a,b^b,c  kA,BkB,c 

\6EIbEIcP22  ^ EIBkBc  +  El  ckB)c  +  EIckA)R 
kAjki 


A,BKB,C 


kA,BkB,c 


Pi  +  1> 


El  c  EIc 

P5  —  1  —  ~j  Pi  5  P6  ~  1  —  7  P2 • 

^5,C  ^5,C 


(A.17) 


(A.18) 


Appendix  B.  Solutions  for  n  equal  length  cracks 


The  solution  for  the  system  of  n  equal  length  cracks  is  determined  by  application  of  the  boundary  con¬ 
ditions  Eq.  (10)  and  continuity  conditions  Eq.  (11).  The  general  solutions  for  the  intact  beam  segment,  with 
index  k  =  0,  and  each  of  the  beam  segments,  with  indices  k  =  1 +  1,  in  the  delaminated  region  using 
the  unconstrained-contact  model  are  given  by  Eq.  (A.2): 

wk{x)  =  CM  +  Ck?x  +  CK  3X2  +  CMx3  (B.l) 


the  4(?z  +  1)  constants  of  integration  are 

c  -c  -0  c  -lPL  c  -  1  p 

Co,,-C0,2-0,  Co, 3--—,  C0A---—, 

P  ( ”+1  1  1  \ 

C,,=.(2  „  + 


Vy_,  £/,  £/„, 

n+1 


Cj 


k,2 


P  n+  \  1  \ 

--(I  +  a)(L-a)  ’ 


y=l  EIj  EI0 / 

n+ 1 


C; 


it, 3 


PL  1  ^  P  l 

c°'4  -  2^  w,'  k~  + 1' 

j= 1  7 


(B.2) 
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The  constants  for  each  of  the  segments  in  the  cracked  region,  C^i,  Q,2,  C^3,  Q?4,  are  identical,  thus  the 
deflections  of  the  beam  segments  in  the  cracked  region  are  the  same.  Therefore  the  other  two  contact 
approximations,  the  constrained-contact  and  the  spring-contact  yield  the  same  solution. 

The  solution  for  the  beam  with  n  equal  length  cracks  could  also  be  determined  by  considering  in  the 
cracked  region  an  equivalent  intact  beam  with  a  reduced  bending  moment  of  inertia  7*  =  •  This  is 

also  true  in  general  under  the  assumptions  of  the  constrained-contact  model  when  the  n  cracks  have  arbi¬ 
trary  lengths.  In  this  case,  the  beam  is  divided  up  into  segments,  each  with  a  reduced  bending  moment  of 
inertia  determined  by  the  sum  of  the  bending  moments  of  inertia  of  the  delaminated  sublaminates  forming 
that  segment. 

The  total  potential  energy  77  of  the  system  is 


1  P2L 3 


6  El  o 


P2a 3 
~6~ 


E 


1 

eTj 


(B.3) 


Appendix  C.  Stability  analysis  of  a  cantilever  beam  with  n  equal  length,  equally  spaced  cracks 


The  stability  of  the  equality  of  length  of  n  equally  spaced  cracks  in  a  cantilever  beam  is  analyzed  here 
using  the  assumptions  of  the  unconstrained-contact  model.  If  a  positive  perturbation  A  a  of  crack  m  is  con¬ 
sidered,  the  energy  release  rate  for  crack  m  is 


ymEh 


6m(N  —  m) 


a  +  A  a\2  /  x3a6  +  6fi  a3  (a  +  A  a)3  +  2x2(a  +  A  a)6^ 

h  )  \  +  x2(a  +  Aaf  , 


(C.l) 


where 


X\  =  6(m4  —  2m3  N  +  m2(N2  —  3)  +  3  mN  —  A2), 
X2  =  3(3 m2  -  3 mN  +  N2), 


*3 


N4 


-  + 


3X1  +  3 Xi  )N 


*i 


N  =  n+  1. 


The  energy  release  rate  for  each  of  the  remaining  n- 1  cracks  is 

fjEh  _  6  TV2  /a\ 2  / y4a6  +  x5a3(a  +  An)3  +  x6(a  +  Aa)6N 
P2  N  -2\h)  y  Xi^3  +  x2(a  +  A  af  , 


(C.2) 


Xa  —  X\  ”*"3^2  "*~^1  + 

Xs  =  - -  Xi )  +  2(z2  +  X1X2  +  3Xi), 
Xe  =  ~\N2{f2  +4X2  +  5zi- 


where 
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The  ratio  of  energy  release  rates,  Eq.  (C.l)  by  Eq.  (C.2),  taking  the  limit  as  Aa  — ►  0  is 

(n2  —  mn  +  2n  —  m  +  1  +  m2)(n  —  1) 


=  m(N  -  m)(N  -  2)  { %3  +  6^  +  3/2 


Nz 


lA  +  15  +  16 


(—m2  +  mn  +  m  —  1  ){n  +  l)2 


(C.3) 


The  equality  of  length  is  stable  if,  for  n  ^  2,  ^m/^i  <  1-  This  condition  is  checked  by  verifying  it  at  the 
extremum  in  terms  of  m  and  at  the  boundaries  for  m=  1,  n.  This  yields 


'  m 


3,  ^ 


<St  n  +  3’ 


1 


(n  +  iy 


(C.5) 


which  are  always  less  than  1  for  n  ^  2. 

If  a  negative  perturbation  Aa  of  crack  m  is  considered  then  the  energy  release  rate  for  crack  m  is 


EimEh  (a  —  Aa\2  1447V 2  a6 

P1  V  *  j  (3 (N  -  2) (a  -  Aa)3  +  (AT  +  6)a3)2  ' 

The  energy  release  rate  for  each  of  the  remaining  n  —  1  cracks  is 

@iEh  g/u\2  XjiiN  +  6)a6  +  6(N  —  2  )a3(a  —  Aa)3]  +  9  %8(a  —  Aa)6 
w  (3(Af  -  2) (a  -  Aa)3  +  (tf  +  6)a3)2 


(C.6) 
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where 


Xl  =  (N2  +  2N  +  3), 

Xs  =  (N3-N  +  2). 

The  ratio  of  energy  release  rates,  Eq.  (C.6)  by  Eq.  (C.7),  taking  the  limit  as  Aa  — ►  0  is 
1441V2  3  n  +  3 

9i  hON+6)+9X%  2n  +  y  ['} 

The  equality  of  length  is  stable  if,  for  n  ^  2,  ^m/^i  >  1,  which  clearly  is  satisfied  by  Eq.  (C.8).  Thus  the 
equality  of  length  is  stable  for  n  equally  spaced  cracks. 


Appendix  D.  Energy  release  rate  expressions  for  a  cantilever  beam  with  two  delaminations 


The  expressions  for  the  normalized  energy  release  rate  for  the  upper  and  lower  delamination,  and 
for  the  cantilever  beam  of  Fig.  5  are  presented  here.  When  the  upper  delamination  is  longer  than  the  lower, 
flu  >  Fig-  5a,  the  solution  shows  that  there  is  contact  along  the  delamination  surfaces.  The  normalized 
energy  release  rates  for  the  upper  and  lower  delamination  have  been  determined  in  closed  form  for  the 
unconstrained-  and  constrained-contact  models.  The  energy  release  rate  for  the  spring-contact  model  is  deter¬ 
mined  through  the  numerical  solution  of  the  problem  and  application  of  Eq.  (9a). 

For  the  constrained-contact  model  the  energy  release  rates  for  a\j  >  a L  are 
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where  H3  =  h^/h  and  H5  =  h5/h.  For  the  unconstrained-contact  model  the  energy  release  rates  for  a\j  >  aL 
are 


Eh  _  fct\j\2  ~  H^)(yx(av/h)6  +  72(au/^)3(^L/^)3  +  73 (^l/^)6) 

-pr~  V  fafa/hf+ysfa/h)3)2 

Eh  =  /aL\2  f(av/h)6H5(l  -  H3f(l  -  H3  -  H5)y6\ 

p 2  V  (yMh)3  +  7s(aL/h)3)2  7’ 

where  the  constants  yi, . .  .  ,y6  depend  on  the  through-thickness  positions  of  the  cracks  and  are 
7!  =  (3 H]  -  3 H3  +  1  ){H\  +  3 H\  +  3HiH5  -  2H3  -  3 Hs  +  l)2, 
y2  =  6H5H\{\  -Hi-  HS){H\  +  3  H\  +  3  H3H5  -  2 H3  -  3  Hs  +  1), 

73  =  3H]H]{\-H3-Hs)2{H]+H3  +  1),  (D.5) 

?4  =  -(3 H\  -  3 Hi  +  1  ){H\  +  3 H]  +  3 H3H5  -  2 H3  -  3 Hs  +  1), 

75  =  -3HsHl(l  -H3-Hs). 

When  the  lower  delamination  is  longer  than  the  upper  delamination,  a^  >  au  Fig.  5b,  the  solution  shows 
that  there  is  opening  along  the  lower  delamination  face,  and  no  contact  along  the  upper  delamination  face, 
namely  the  two  upper  beam  segments  have  the  same  vertical  deflection.  Therefore,  the  spring-  and  uncon¬ 
strained-contact  models  lead  to  the  same  solution.  The  constrained-contact  model,  which  prevents  opening, 
leads  to  a  solution  that  is  clearly  incorrect  in  this  regime.  However,  as  will  be  seen  later,  it  is  a  bound  of  the 
exact  solution,  which  is  useful  for  delamination  configurations  with  more  delaminations  where  the  exact 
solution  is  not  so  obvious.  The  energy  release  rates  are  given  by  Eqs.  (D.1)-(D.5)  with  the  indices  of  crack 
positions  and  lengths  reversed,  aL  — >  a\j ,  au  — »  aL,  — >  H5  and  H5  — >  H3. 

When  the  delaminations  have  the  same  length,  the  energy  release  rate  of  each  delamination  when  the 
delaminations  propagate  simultaneously,  ^B,  is  given  by  Eq.  (12). 

As  for  the  case  of  equal  length  delaminations,  Fig.  3,  maps  can  be  constructed  for  a  general  through¬ 
thickness  distribution  and  lengths  of  the  delaminations  that  show  regions  in  which  the  energy  release  rate 


(D.3) 

(D.4) 


Fig.  D.l.  Maps  of  regions  of  different  energy  release  rates  for  a  system  of  two  unequal  length  cracks  (schematic  of  Fig.  5): 
(a)  unconstrained  model,  (b)  constrained  model. 
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of  one  of  the  cracks  is  higher  or  lower  than  that  of  the  other  crack.  Maps  for  the  constrained-  and  uncon¬ 
strained-contact  models  are  shown  in  Fig.  D.l. 

The  maps  depend  on  the  through-thickness  position  of  the  two  cracks  and  incorporate  contours  of  equal 
^  of  both  cracks.  The  contours  depend  on  the  relative  length  of  the  two  cracks.  Dashed  contours  corre¬ 
spond  to  Ajj  =  (ajj  —  <2l)A*u  and  are  to  be  used  when  au  >  aL.  Solid  contours  correspond  to 
AL  =  (<2l  -  au)/a^  and  are  to  be  used  when  aL  >  ajj.  The  map  in  Fig.  D.la  has  been  constructed  using 
the  unconstrained-contact  model  and  the  map  in  Fig.  D.lb  using  the  constrained- contact  model.  A  more 
complicated  map,  depending  also  on  the  length  of  the  cracks,  can  be  constructed  using  the  spring-contact 
model. 

The  shaded  region  in  each  figure  refers  to  cracks  of  the  same  length,  Av  =  AL  =  0,  and  define  configu¬ 
rations  for  which  the  energy  release  rate  of  each  crack  when  the  cracks  propagate  simultaneously,  given 
by  Eq.  (12),  is  maximum  (the  shaded  region  in  Fig.  D.la  coincides  with  that  of  Fig.  3).  Points  to  the  left  or 
above  the  contour  corresponding  to  the  relative  length  of  the  two  cracks  define  through-thickness  distribu¬ 
tions  of  the  cracks  for  which  the  upper  crack  has  the  higher  energy  release  rate.  Points  below  or  to  the  right 
of  the  contour  define  through- thickness  distributions  of  the  cracks  for  which  the  lower  crack  has  the  higher 
energy  release  rate.  The  evolution  of  the  cracks  can  be  followed  in  this  diagram  by  updating  the  contour 
from  that  corresponding  to  the  current  A,  to  the  configuration  reached  by  the  last  crack  growth  event. 
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ABSTRACT 

The  problem  is  considered  of  a  fiber  that  is  driven  dynamically,  by  compression  at  one 
end,  into  a  matrix.  The  fiber  is  not  initially  bonded  to  the  matrix,  so  that  its  motion  is 
resisted  solely  by  friction.  Prior  work  based  on  simplified  models  has  shown  that  the 
combination  of  inertial  effects  and  friction  acting  over  long  domains  of  the  fiber-matrix 
interface  give  rise  to  behaviour  that  is  far  more  complex  than  in  the  well-known  static 
loading  problem.  The  front  velocity  may  depart  significantly  from  the  bar  wave  speed 
and  regimes  of  slip,  slip/stick,  and  reverse  slip  can  exist  for  different  material  choices  and 
loading  rates.  Here  more  realistic  numerical  simulations  and  detailed  observations  of 
dynamic  displacement  fields  in  a  model  push-in  experiment  are  used  to  seek  more 
complete  understanding  of  the  problem.  The  prior  results  are  at  least  partly  confirmed, 
especially  the  ability  of  simple  shear  lag  theory  to  predict  front  velocities  and  gross 
features  of  the  deformation.  Some  other  fundamental  aspects  are  newly  revealed, 
including  oscillations  in  the  interface  stresses  during  loading;  and  suggestions  of 
unstable,  possibly  chaotic  interface  conditions  during  unloading.  Consideration  of  the 
experiments  and  two  different  orders  of  model  suggest  that  the  tentatively  characterized 
chaotic  phenomena  may  arise  because  of  the  essential  nonlinearity  of  friction,  that  the 
shear  traction  changes  discontinuously  with  the  sense  of  the  motion,  rather  than  being 
associated  with  the  details  of  the  constitutive  law  that  is  assumed  for  the  friction.  This 
contrasts  with  recent  understanding  of  instability  and  ill-posedness  at  interfaces  loaded 
uniformly  in  time,  where  the  nature  of  the  assumed  friction  law  dominates  the  outcome. 
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1. 


Introduction 


This  paper  will  address  certain  aspects  of  the  more  general  problem  illustrated  in  Figure 
1.  In  the  general  problem,  a  fiber  that  is  initially  bonded  to  a  matrix  is  loaded 
dynamically  in  either  tension  or  compression  at  its  end,  with  a  load  that  is  some  function 
of  time.  The  fiber  and  matrix  are  elastically  dissimilar.  The  assembly  may  be  of  plane  or 
axisymmetric  geometry.  Boundary  conditions  for  the  axisymmetric  case  may  be  those  of 
Type  I  or  Type  II  (i.e.,  stress  or  displacement  conditions  at  the  boundary  of  a  cylindrical 
cell),  as  designated  by  Hutchinson  and  Jensen  (Hutchinson  and  Jensen,  1990)  to  simulate 
an  array  of  fibers;  or  those  corresponding  to  an  isolated  fibre  in  a  semi-infinite  body. 
Boundary  conditions  for  a  plane  problem  may  be  periodic  to  simulate  an  array;  or  again 
those  of  an  isolated  fiber  in  a  semi-infinite  body.  A  state  of  initial  compression  may  exist 
across  the  fiber-matrix  interface,  e.g.,  due  to  residual  thermal  stresses.  The  load  causes 
the  fiber  to  debond  progressively  from  the  matrix,  allowing  relative  sliding  between  the 
two,  which  is  opposed  by  friction.  The  constitutive  law  of  the  friction  may  involve 
displacement,  displacement  rate,  and  the  magnitude  of  the  interface  compression.  The 
debonding  mechanism  may  result  in  a  nonlinear  process  zone  in  the  debond  crack  wake 
of  significant  length.  Friction  will  generally  act  over  much  longer  lengths,  often  the 
whole  domain  of  sliding. 

In  this  paper,  model  experiments  and  numerical  simulations  of  the  push-in  problem 
(compressive  end  loading)  will  be  used  to  study  the  particular  limit  that  the  interface  is 
initially  debonded  (limit  of  zero  debond  toughness).  One  motivation  is  to  explore  the 
connections  between  the  physics  of  long  range  frictional  sliding  and  shear  fracture,  using 
model  experiments  in  which  a  reduced  set  of  mechanisms  operates.  Simultaneously,  the 
question  is  addressed  of  how  simply  dynamic  interface  failure  in  composites  can  be 
represented,  without  loss  of  accuracy.  The  latter  query  is  inspired  by  watershed 
experiments  in  the  static  case  on  cylindrical  fibers  in  a  matrix  (e.g.,  (Marshall,  1992)).  In 
the  static  case,  shear  lag  analysis  of  the  fiber/matrix  debonding  problem,  which  assumes 
Lame-like  field  solutions,  proves  very  accurate  for  most  parameter  ranges  (Hutchinson 
and  Jensen,  1990)  and  forms  an  immediate  link  to  the  problem  of  composite  fracture 
(Budiansky  et  al.,  1986;  Marshall  et  al.,  1985).  Where  static  shear  lag  models  are 
accurate,  analytical  results  usually  stand  independently  of  whether  the  fibers  are  more  or 
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less  compliant  than  the  matrix  and  whether  plane  or  axisymmetric  conditions  exist; 
preliminary  shear  lag  analyses  of  the  dynamic  case  suggest  the  same  compass,  although 
the  physics  are  considerably  more  complicated  than  in  the  static  case  (Sridhar  et  al., 
2003).  The  problem  of  Fig.  1  may  also  be  considered  as  a  particular  case  of  delamination 
failure  in  a  symmetrically  laid  up  laminated  composite,  for  which  the  plane  geometry  is 
directly  applicable  and  the  lamina  represented  by  the  fiber  may  be  more  or  less  compliant 
than  the  other  layers.  The  questions  of  how  accurate  shear  lag  theory  might  be  in 
dynamic  cases,  and  what  subset  of  the  physics  of  interface  failure  might  remain 
reasonably  well  represented  in  it,  are  therefore  of  practical  interest.  Issues  specific  to  the 
push-in  problem  will  be  addressed  here,  including  how  the  degrees  of  freedom  used  in 
representing  displacement  fields  in  modeling  and  the  assumed  nature  of  the  friction  law 
influence  stability,  front  formation,  front  velocity,  stress  fields,  etc. 

Friction  Acting  Along  Bars 

Shear-lag  models  of  the  dynamic  case  of  the  problem  of  Fig.  1,  with  friction  that  is 
spatially  and  temporally  uniform  in  magnitude  (but  not  in  sign),  have  revealed  the 
regimes  of  behaviour  that  might  be  expected  for  an  initially  debonded  interface  (Cox  et 
al.,  2001;  Nikitin  and  Tyurekhodgaev,  1990;  Sridhar  et  al.,  2003).  First,  the  front  of  the 
furthest  deformation  in  the  fiber  does  not  travel  at  the  longitudinal  wave  speed  for  the 
fiber,  Cf,  but  at  some  other  velocity  that  depends  on  the  friction  strength,  the  elastic 
mismatch,  the  wave  speeds  in  the  fiber  and  the  matrix,  and  the  loading  rate.  The  front 
velocity  is  bounded  by  Cf  and  the  longitudinal  wave  speed  in  the  matrix,  Cm,  and  may  be 
lower  or  higher  than  Cf.  Second,  depending  on  the  same  parameters,  distinct  regimes  of 
slip,  slip  and  stick,  and  reverse  slip  exist,  even  when  the  load  point  displacement  is 
monotonic.  Reverse  slip  refers  to  the  condition  that  the  relative  motion  of  the  fiber  and 
matrix  is  opposite  in  sense  to  that  of  the  load  point.  Reverse  slip  is  possible  when  the 
wave  speed  in  the  matrix  is  higher  than  that  in  the  fiber,  so  that  the  fiber  is  pulled  along, 
at  the  deformation  front,  by  the  matrix.  The  domains  of  slip,  stick,  and  reverse  slip  are 
fixed  (growing  similarly  with  time)  for  linearly  increasing  loads,  but  appear  and 
disappear  in  complicated  sequences  for  other  loading  cases,  even  those  that  are  simple 
functions  of  time  (e.g.,  step  functions).  The  solutions  are  much  more  complicated  than 
for  the  same  problem  in  static  loading,  where  the  only  history  dependence  that  affects  the 
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solution  is  the  order  of  any  load  reversals  [e.g.,  (Marshall,  1992)].  There  is  even  a  hint 
(unproven)  that  the  birth  and  death  of  slip,  reverse  slip,  and  stick-slip  domains  may  be 
chaotic,  with  arbitrarily  small  changes  in  the  load  point  history  resulting,  at  later  times,  in 
finitely  dissimilar  domain  patterns.  Numerical  work  has  supported  the  shear  lag  results 
for  linearly  increasing  loading,  under  the  modeling  conditions  assumed,  for  most 
parameter  values  (Sridhar  et  al.,  2003).  However,  the  strong  non-linearity  of  the  friction 
law,  which  may  reverse  sense  at  moving  boundaries,  makes  accurate  numerical  location 
of  the  boundaries  very  challenging,  even  under  the  constraint  of  assumed  Lame-like 
solutions  (Sridhar,  Dunn,  and  Cox,  unpublished  work,  2002). 

Friction  Between  Half-Spaces 

The  problem  of  an  interface  with  no  debond  energy  but  long,  possibly  infinite,  domains 
of  friction  has  been  studied  extensively  in  the  geometry  of  two  remotely  and  uniformly 
loaded  half-spaces,  contacting  on  a  plane.  Among  a  number  of  insightful  articles,  the 
introductory  sections  of  (Cochard  and  Rice,  2000)  provide  an  especially  helpful 
summary.  Some  salient  results  are  as  follows.  1)  For  elastically  dissimilar  materials 
subject  to  Coulomb  friction,  uniform  slip  becomes  unstable  and  slip  pulses  can  form 
instead,  i.e.,  slip  domains  on  either  side  of  which  the  interface  is  not  sliding  (Adams, 
1998;  Weertman,  1980).  Slip  pulses  can  propagate  in  either  direction,  in  different  cases, 
and  at  values  of  the  remote  shear  stress  that  are  less  than  the  friction  stress.  2)  Unstable 
slip  is  also  possible  for  elastically  homogeneous  cases,  but  not  for  a  simple  Coulomb  law. 
3)  Cochard  and  Rice  go  on  to  distinguish  instability  from  ill-posedness:  the  former  exists 
when  spatial  perturbations  to  the  state  of  slip  grow  in  amplitude;  the  latter  refers  to  the 
non-existence  of  solutions  of  any  kind.  Instability  does  not  preclude  the  possibility  of 
uniform  slip  as  a  formal  solution,  given  the  hypothetical  presence  of  perfectly  uniform 
conditions.  Ill-posedness  points  to  an  inconsistency  in  the  physics  of  the  problem  as 
stated.  4)  Whether  the  interface  problem  exhibits  ill-posedness  or  instability  depends  on 
both  the  degree  of  elastic  mismatch  and  the  nature  of  the  assumed  friction  law.  The 
relevant  measure  of  elastic  mismatch  is  the  point  at  which  the  generalized  Rayleigh 
velocity  (a  mode  of  interface  wave  propagation  found  for  a  frictionless  interface  along 
which  loss  of  contact  does  not  occur)  ceases  to  be  defined.  Ill-posedness  prevails  for  all 
friction  coefficients  when  the  elastic  mismatch  is  mild  (generalized  Rayleigh  wave 
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exists);  and  for  friction  coefficients  above  a  modest  critical  value  when  the  mismatch  is 
more  severe.  Thus  ill-posedness  is  the  common  case.  5)  Ranjith  and  Rice  (Ranjith  and 
Rice,  2001)  assign  blame  to  the  simplified  physics  of  Coulomb’s  law,  in  particular  its 
implication  that  friction  can  change  discontinuously  in  time  upon  a  finite  change  in  the 
normal  contact  pressure.  A  modified  law,  in  which  the  change  occurs  smoothly  over  a 
characteristic  response  time,  had  in  fact  already  been  inferred  from  experiments  (Prakash 
and  Clifton,  1993;  Prakash,  1998).  Such  a  law  regularizes  the  ill-posed  planar  interface 
problem;  and  the  regularized  solutions  continue  to  exhibit  slip-pulse  character  (Cochard 
and  Rice,  2000). 

The  complexity  seen  in  the  Lame-like  solutions  to  the  push- in  problem  of  Fig.  1  has 
different  physical  origins  to  the  instability  and  ill-posedness  discussed  in  the  last  two 
paragraphs.  In  the  work  based  on  Lame-like  solutions,  the  friction  force  was  assumed  to 
be  uniform  and  constant,  apart  from  sign  reversals,  and  not  influenced  by  changes  in  the 
normal  compression  at  the  interface.  Complexity  arises  from  the  effects  of  friction  on  the 
propagation  of  wave  pulses  along  the  fiber  and  matrix  system,  excited  by  the  dynamic 
end  loading,  and  is  associated  with  the  birth  and  death  of  slip,  stick,  and  reverse  slip 
domains,  which  is  complicated  by  the  strength  of  the  nonlinear  effects  due  to  the  sign- 
reversal  property  of  the  friction.  In  the  problem  of  planar  frictional  interface  between 
half-spaces,  instability  and  ill-posedness  are  direct  consequences  of  the  assumed  relation 
between  friction  and  the  normal  traction,  and  depend  strongly  on  whether  Coulomb’s  or 
another  law  is  assumed. 

The  final  point  from  the  body  of  work  on  planar  interfaces  that  is  of  present  relevance  is 
the  relation  of  the  boundedness  of  slip-pulse  domains  to  the  nature  of  the  assumed 
friction  law.  Zheng  and  Rice  (Zheng  and  Rice,  1998)  looked  at  this  question  with  a 
velocity-weakening  friction  law,  expressed  as  a  decrease  in  the  coefficient  of  friction  as 
the  shear  displacement  rate  rises.  (This  law  is  similar  to  that  of  Prakash  and  Clifton,  but 
expressed  in  inverse  form,  in  terms  of  the  displacement  rate  rather  than  the  rate  of  change 
of  the  shear  traction.)  For  such  a  law,  slip  motion  can  occur  either  as  a  crack-like  event, 
in  which  the  interface  slips  simultaneously  over  a  semi-infinite  domain  extending  back 
from  a  slip  front;  or  as  a  slip-pulse  motion,  with  slip  confined  to  finite,  moving  domains. 
Which  case  prevails  (in  the  problem  of  contacting  half-spaces)  depends  on  the  rate  of 


4 


velocity  weakening:  if  the  slope  of  the  decline  is  small,  crack-like  slip  behaviour  occurs; 
if  it  is  large,  slip-pulse  motion  occurs  (Zheng  and  Rice,  1998). 

Effects  of  Interface  Bonding 

A  separate,  major  body  of  literature  has  been  written  on  the  problem  of  interface  shear 
crack  propagation  in  systems  in  which  the  interface  is  initially  bonded.  The  most 
pertinent  results  here  refer  to  the  question  of  allowed  and  preferred  crack  velocities. 
Energy  considerations  suggest  that,  unlike  mode  I  cracks,  which  are  forbidden  above  the 
shear  wave  speed,  shear  cracks  can  in  principle  propagate  at  any  velocity,  c,  in  the  so- 
called  intersonic  regime,  Cs  <  c  <  Cl,  where  Cs  and  Cl  are  the  shear  and  longitudinal 
wave  speeds  (Broberg,  1994;  Burridge  et  al.,  1979).  However,  numerical  studies  imply  a 
restrictive  condition  for  intersonic  propagation,  that  the  crack  tip  process  zone  be  diffused 
over  a  finite  interval  and  not  concentrated  at  a  point,  as  in  a  classical  brittle  crack 
(Andrews,  1976).  In  the  work  of  Burridge  et  al.,  the  velocity  V2CS  assumes  a  special 
role:  for  c  <  Cs,  crack  acceleration  is  unstable  (associated  with  decreasing  load);  while 
for  c  >  V2Cs,  crack  acceleration  is  stable  (increasing  load). 

These  early  theoretical  deductions  have  now  been  complemented  by  experiments  that 
have  achieved  shear  crack  conditions  for  the  first  time  in  the  laboratory  by  exploiting  the 
propensity  of  laminated  materials  to  confine  cracks  to  prescribed  planes  (Rosakis  et  al., 
1999).  Specimens  consisting  of  monolithic  slabs  laminated  at  a  bonded  interface  as  well 
as  composite  ply  laminates  have  been  studied.  Under  appropriate  dynamic  loading, 
approximately  mode  II  conditions  can  be  maintained,  in  contrast  to  tests  with  a  uniform 
body  (free  of  pre-existing  weak  planes),  where  crack  deflection  or  bifurcation  would 
quickly  cause  a  transition  to  mode  I  conditions.  The  laminate  experiments  confirm  mode 
II  crack  propagation  in  the  intersonic  regime  and  are  also  consistent  with  the  predicted 
stability  transition  at  V2Cs,  since  the  crack  is  observed  to  accelerate  rapidly  past  the  shear 
wave  speed,  but  decelerate  and  achieve  approximately  uniform  velocity  at  V2Cs. 
Needleman  has  performed  numerical  simulations  of  this  particular  experimental 
configuration,  where  the  driving  force  is  a  dynamic  impact  end-load  rather  than  the 
remote  shear  loading  of  the  early  theoretical  studies,  and  confirmed  most  aspects  of  the 
observed  behaviour  (Needleman,  1999).  Needleman’s  model  incorporates  mixed  mode 
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cohesive  elements  to  represent  the  debonding  process,  with  cohesive  laws  (traction  vs. 
displacement)  that  feature  both  an  initial  hardening  phase  and  a  decaying  tail.  These  laws 
are  somewhat  different  from  that  in  Burridge  et  al.,  which  was  a  purely  softening  law  and 
for  mode  II  displacements  only.  The  details  of  the  law  do  affect  the  outcome  in  some 
ways.  For  example,  while  Needleman’s  simulations  show  the  same  qualitative  history  of 
rapid  acceleration  to  a  nearly  constant  velocity  as  seen  in  experiments  (Rosakis  et  al., 
1999),  the  value  of  the  final  velocity  generally  exceeds  V2Cs,  rather  than  equalling  it, 
with  the  speed  attained  depending  on  the  maximum  stress  assumed  for  the  cohesive  law. 
In  a  separate  theoretical  study  of  shear  crack  propagation  along  bi-material 
(Homalite/steel)  interfaces,  a  large  jump  in  terminal  crack  velocity  was  found  at  a  certain 
prescribed  load  point  velocity  (Needleman  and  Rosakis,  1999).  For  inferior  load  point 
velocities,  crack  propagation  was  limited  by  the  Rayleigh  wave  speed  in  the  Homalite 
(the  slower  medium);  above  the  transition,  the  terminal  velocity  approximately  doubled, 
to  be  near  the  longitudinal  wave  speed  of  the  Homalite. 

While  intersonic  propagation  is  physically  admissible,  shear  crack  propagation  remains 
forbidden  in  the  velocity  interval,  Cr  <  c  <  Cs,  where  Cr  is  the  Rayleigh  velocity.  The 
mechanism  by  which  a  crack  can  jump  across  this  velocity  gap  involves  the  creation  of  a 
daughter  crack  ahead  of  the  parent  crack,  which  then  links  to  the  parent  in  an  effective 
surge  of  crack  growth  (Andrews,  1976;  Burridge,  1973;  Gao  et  al.,  2001).  While 
relationship  of  this  phenomenon  to  slip-pulse  solutions  for  planar  frictional  interfaces 
must  exist,  the  presence  of  a  velocity  gap  must  be  associated  with  presence  of 
concentrations  of  tractions  where  interface  debonding  is  taking  place,  which  do  not 
generally  appear  in  interfaces  coupled  solely  by  friction. 

Contrasts  in  Physics 

Thus  shear  crack  propagation  analysis,  pulse-slip  solutions  for  the  planar  frictional 
interface,  and  Lame-like  solutions  of  the  push-in  problem  with  zero  debond  energy  (Fig. 
1)  all  show  that  variable  front  velocities  up  to  the  maximum  of  the  longitudinal  wave 
speeds  in  the  fiber  and  the  matrix  are  physically  admissible.  However,  important 
distinctions  arise  between  the  different  problems.  In  particular,  there  is  no  analogue  in 
the  Lame-like  solutions  of  the  push-in  problem  with  zero  debond  energy  to  the  velocity 
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gap  in  the  interval  (Cr,  Cs)  in  fracture  problems,  nor  any  special  role  for  the  velocity 
V2Cs.  Neither  do  direct  analogues  of  the  phenomena  of  instability  and  ill-posedness  that 
have  been  so  extensively  studied  for  planar  frictional  interfaces  appear  in  the  existing 
Lame-like  solutions,  since  these  have  been  developed  for  uniform  and  constant  (not 
Coulomb)  friction. 

Solutions  in  the  literature  for  dynamic  shear  cracks  and  the  initially-debonded  push-in 
problem  also  differ  in  another  very  important  way.  Even  though  cohesive  models  have 
been  used  for  dynamic  shear  crack  modeling,  thus  spreading  the  debond  process  into  a 
diffuse  process  zone,  the  dimensions  of  this  process  or  cohesive  zone  have  been  relatively 
small.  For  example,  in  the  simulations  of  Needleman  (Needleman,  1999),  the  cohesive 
zone  is  ~  1  mm  (a  result  of  the  assumed  parameter  values  in  the  cohesive  law),  while  the 
crack  propagates  over  distances  of  ~  25  mm.  Behind  the  cohesive  zone,  the  fracture 
surfaces  remain  traction  free.  A  small  process  zone  representing  debonding  is  very 
different  from  the  conditions  expected  in  the  presence  of  friction,  especially  in  mode  II 
problems,  where  shear  tractions  may  act  over  the  whole  cracked  specimen.  While  a  shear 
cohesive  traction  enters  the  fracture  problem  in  exactly  the  same  way  that  a  frictional 
traction  would  enter,  the  extension  of  cohesive  modeling  to  friction  zones  that  may  be  as 
long  as  the  crack  has  not  been  made  in  dynamic  fracture  studies.  (Calculations  for  bi¬ 
material  interfaces  did  predict  the  existence  of  large-scale  contact  zones  trailing  the 
crack-tip  (Needleman  and  Rosakis,  1999),  but  friction  in  these  zones  was  not  modelled.) 
Models  have  either  treated  long  zones  of  friction  in  the  absence  of  a  debond  process  (the 
Lame-like  push-in  problem  or  the  planar  frictional  interface  problem);  or  debonding  as  a 
relatively  short  process  zone  in  the  absence  of  long  zones  of  friction. 

Final  introductory  comments  are  directed  upon  the  nature  of  the  cohesive  law.  In  static 
fracture  of  a  material  that  is  not  rate  dependent,  only  two  or  three  degrees  of  freedom  in 
the  cohesive  law  have  a  significant  influence  on  the  fracture  response  (e.g.,  (Bao  and  Suo, 
1992;  Cox  and  Marshall,  1994;  Massabo  and  Cox,  1999)).  If  the  zone  has  finite  extent, 
most  characteristics  of  crack  propagation  are  determined  by  two  parameters,  which  may 
be  taken  as  the  critical  traction  in  the  cohesive  law  and  the  area  under  the  law  (work  of 
fracture).  If  the  cohesive  law  incorporates  an  initial  hardening  phase  (traction  rising  with 
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displacement),  then  the  cohesive  mechanism  may  act  over  the  whole  crack,  even  for  very 
long  cracks,  in  which  case  a  third  parameter,  the  slope  of  the  hardening  phase,  controls 
crack  growth  (Cox  and  Marshall,  1994).  If  the  cohesive  mechanism  exhibits  rate 
dependence  (excluding  inertial  effects),  e.g.,  due  to  creep  or  viscosity,  then  one  or  at  most 
two  additional  parameters,  describing  the  rate  of  decay  of  the  cohesive  tractions,  appear 
to  suffice  to  predict  fracture  behaviour  (Cox  and  Sridhar,  2001).  The  possibilities  when 
inertial  effects  enter  may  be  richer.  For  example,  in  none  of  the  work  on  static  fracture 
are  there  cases  where  the  rate  of  decline  of  the  law  with  increasing  displacement  plays  a 
significant  role  (except  as  it  affects  the  area  under  the  curve,  i.e.,  the  work  of  fracture). 
Yet  in  dynamic  studies  of  friction,  the  rate  of  decline  of  the  coefficient  of  friction,  i.e.,  of 
the  cohesive  shear  traction  due  to  friction,  with  displacement  rate  has  a  profound  effect:  it 
determines  whether  slip  occurs  in  a  crack-like  or  slip-pulse  mode  (Zheng  and  Rice, 
1998).  This  and  other  results,  including  theoretical  models  (Lapusta  and  Rice,  2003)  and 
novel  laboratory  simulations  of  slip-pulse  motion  (Xia  et  al.,  2004),  show  that 
constitutive  laws  for  dynamic  friction  need  to  be  represented  with  more  degrees  of 
freedom  than  laws  in  static  facture.  Inertial  effects  appear  to  manifest  further  details  of 
the  traction  law  in  fracture  behaviour. 

All  of  these  background  works  have  addressed  some  aspects  of  the  general  problem  of 
Fig.  1,  but  they  remain  separate  pieces  of  understanding,  not  yet  connected  to  one 
another.  Thus  large-domain  friction  problems  have  not  yet  embraced  the  influence  of  a 
strong  bond  (energetically  significant  crack  tip  process);  while  fracture  work  has  not  yet 
extended  to  very  large  (possibly  semi-infinite)  domains  of  friction.  Generality  in  the 
friction  law  has  begun  to  appear  only  in  studies  devoted  to  friction  as  a  phenomenon,  but 
these  studies  imply  that  over-simplifying  friction  laws  can  have  a  strong  affect  on  the 
possible  behaviour  in  other  systems.  And  last,  the  most  general  studies  of  fiber  push-in 
as  a  special  problem,  with  general  end-loading  histories,  remain  to  be  undertaken. 

2.  Experiments 

A  model  planar  specimen  geometry  was  designed  to  investigate  dynamic  fiber  push-in  for 
a  bi-material  system  in  which  no  chemical  bond  exists  at  the  material  interface.  A 
transparent  and  birefringent  polymer  sheet  (Homalite-100),  acting  as  a  “fiber,”  was 
placed  between  two  steel  plates,  analogous  to  a  “matrix”  (Fig.  2).  The  dynamic  fiber 
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push- in  process  was  mimicked  by  loading  the  fiber  with  a  projectile  traveling  at  speeds 
ranging  from  10  to  40  m/sec.  The  resistance  of  the  polymer/metal  interfaces  to  slip  was 
manipulated  by  applying  a  static  compressive  pre-load  acting  perpendicular  to  the 
interfaces,  whose  value  varied  from  0.5  -  175  MPa. 

Discussion  of  the  resulting  deformation  will  refer  to  the  coordinate  system,  (xi,  xi),  with 
the  origin  located  at  the  mid-plane  of  the  Homalite  piece,  at  the  end  where  it  is  impacted 
(Figs.  1  and  2).  The  half-height  of  the  Homalite  will  be  denoted  h.  Thus  the  interfaces 
lie  at  X2  =  ±h.  The  origin  of  time,  t,  is  the  moment  of  first  impact. 

Materials  and  Specimens:  The  densities  and  relevant  wave  speeds  of  the  test  materials 
are  given  in  Table  1.  All  plates  were  9.5  mm  thick  and  82.5  mm  in  length.  The  steel 
pieces  were  50.8  mm  in  height.  The  height,  2 h,  of  the  Homalite  piece  was  16.5  mm.  The 
Homalite  dimensions  were  chosen  to  provide  a  fiber  aspect  ratio,  defined  as  length 
divided  by  the  half  height,  h,  of  10.  The  steel  height  was  chosen  to  minimize  the 
interactions  of  wave  reflections  from  the  specimen  boundaries  with  the  dynamic  sliding 
process.  This  assures  that  the  case  studied  is  equivalent  to  that  of  an  isolated  fiber  in  a 
semi-infinite  body.  The  pieces  of  the  specimen  were  aligned  and  stacked  vertically 
(steel/Homalite/steel)  without  the  use  of  adhesive  or  bonding  agent. 


Table  1.  Density  and  elastic  wave  speeds  for  the  test  materials 


density,  p 
(kg/m3) 

Cr 

(m/s) 

Cs 

(m/s) 

Cl 

(m/s) 

Homalite 

1230 

1087 

1187 

2060 

steel 

8000 

2977 

3254 

5443 

Dynamic  Loading  and  Characterization:  The  fiber  was  loaded  dynamically  by  impacting 
its  end  with  a  25.4  mm  diameter  cylindrical  steel  projectile,  44.5  mm  in  length,  which 
was  accelerated  using  a  light  air  gun.  Projectile  speeds  were  varied  over  the  range  10  - 
40  m/s  by  varying  the  gun  pressure.  A  small  steel  tab  was  affixed  to  one  end  of  the 
Homalite  piece  to  eliminate  the  possibility  of  impact  damage  to  the  end  of  the  relatively 
brittle  polymer  during  dynamic  loading  (Fig.  2).  The  dynamic  sliding  of  the  Homalite 
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relative  to  the  steel  was  characterized  using  high-speed  photography  in  conjunction  with 
dynamic  photoelasticity.  A  strain  gage  was  mounted  to  the  surface  of  the  steel  tab  to 
provide  a  trigger  signal  for  the  initiation  of  the  recording  of  high-speed  photographs 
(described  below).  From  sequences  of  images,  the  dynamic  loading  history,  o(t),  of  the 
fiber,  the  initiation  time  for  the  onset  of  interface  sliding,  and  the  interface  slip  speed 
histories  were  obtained. 

Dynamic  photoelasticity  apparatus  was  configured  as  follows.  A  collimated  laser  beam, 
100  mm  in  diameter,  was  passed  through  a  circular  polarizer,  the  length  of  the  Homalite, 
and  a  second  polarizer  to  generate  the  photoelastic  images,  which  were  recorded  using  a 
high-speed  digital  camera  (Cordin  220-16,  Cordin  Scientific  Imaging,  Salt  Lake  City, 
Utah,  84119).  The  camera  was  capable  of  recording  16  frames.  The  frames  were 
recorded  every  few  microseconds;  the  interframe  times  were  varied  in  each  experiment. 
In  dynamic  photoelasticity,  the  fringes  show  constant-value  contours  of  the  difference  of 
principal  stresses,  01-02. 

2.1  Results 

General  observations: 

A  total  of  16  experiments  were  conducted.  The  role  of  impact  speed  (loading  rate)  was 
investigated  at  a  constant  static  pre-load;  and  that  of  pre-load  was  investigated  at  a 
constant  impact  speed.  Figure  3  shows  a  representative  high-speed  isochromatic  fringe 
pattern  recorded  during  sliding  in  the  common  member  of  these  two  test  sets  ( cr®  =  50 

MPa,  vo  =  38±2  m/s).  The  approximately  vertical  fringes  represent  the  stress  generated 
from  the  impact  (from  the  left)  and  propagate  to  the  right.  Shortly  behind  the  propagating 
front,  kinks  form  in  the  fringes  in  angled  bands.  These  bands  intersect  the  top  and  bottom 
interfaces,  implying  non-smooth  variations  in  the  interface  displacement.  Among  the 
many  fringe  patterns  collected  from  all  tests,  consistency  in  the  pattern  of  bands  was 
imperfect.  In  some  images,  either  noise  or  experimental  variance  made  it  difficult  to 
confirm  or  refute  the  presence  of  bands  that  might  have  been  expected  to  be  present 
based  on  their  presence  in  other  images.  In  the  images  with  better-defined  fringes,  two 
bands  were  prominent,  typified  by  the  dashed  white  lines  in  Fig.  3.  Based  on  changes  in 
the  spacing  or  angle  of  fringes  before  and  after  a  band,  or  in  the  derivative  of  the  spacing, 
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the  first  of  these  is  interpreted  as  the  beginning  of  a  domain  of  increasing  interface  shear 
stress;  the  second  as  the  point  of  transition  to  interfacial  sliding.  In  Fig.  3,  the  domain 
marked  “increasing  friction”  (or  interfacial  shear  stress)  is  clearly  defined  as  such  by  the 
fringe  behavior;  in  the  domain  marked  “constant  friction,”  the  fringe  variations  are 
consistent  with  this  interpretation,  but  the  implication  is  not  as  clear  due  to  noise.  A  peak 
in  the  fringes  on  the  specimen  mid-plane  locates  the  onset  of  unloading. 

Figure  4  presents  images  similar  to  that  in  Fig.  3  but  taken  from  experiments  conducted 
with  higher  and  lower  static  pre-loads.  The  difference  in  the  visibility  of  bands  in  these 
figures  and  Fig.  3  is  representative.  Further  inspection  of  Figs.  3  and  4  shows  that  as  the 
pre-load  increases,  the  kinking  of  the  fringes  is  associated  with  larger  distortions, 
suggesting  that  the  onset  of  sliding  is  associated  with  a  sharper  discontinuity.  A  similar 
trend  was  observed  as  the  impact  speed  was  increased.  The  dynamic  sliding  process  is 
associated  with  higher,  more  localized  stresses  at  higher  loading  rates  and  higher  static 
pre-load. 

Sliding  Speeds: 

The  velocity  of  the  dominant  band  in  each  image  along  the  interface  was  estimated  by 
comparing  its  position  in  successive  frames.  Figure  5  shows  plots  of  the  positions  of  the 
intersections  of  a  symmetric  pair  of  such  bands  with  the  top  and  bottom  interfaces  as  a 
function  of  time  for  an  experiment  conducted  with  12.5  MPa  pre-load  and  38  m/s  impact 
speed.  The  bands  propagate  at  the  same  velocity  on  the  top  and  bottom  and  the  variation 
of  sliding  position  with  time  is  essentially  linear.  This  contrasts  with  the  common  case 
for  shear  crack  propagation,  where  a  period  of  acceleration  of  the  crack  tip  is  often 
observed,  which  may  approach  constant  velocity  asymptotically.  Here,  the  velocity  of 
the  bands  is  constant  throughout  the  experiment,  within  the  experimental  error  of 
approximately  ±3%  (typically  on  the  order  of  ±50  m/s).  These  observations  were 
consistent  in  all  of  the  experiments,  indicating  that  on  this  scale  of  observation,  sliding  is 
a  steady-state  process.  In  this  particular  experiment,  the  slope  of  the  plot  in  Fig.  5 
indicates  a  sliding  speed  of  2050  m/s  or  approximately  the  dilatational  wave  speed  of 
Homalite.  Model  analysis  (see  below)  suggests  that  these  bands  probably  corresponded 
to  the  beginning  of  a  domain  of  enhanced  interfacial  shear  stress. 
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Figure  6  provides  a  summary  of  the  velocities  of  the  dominant  bands  measured  from  all 
of  the  experiments,  plotted  as  a  function  of  static  pre-load.  While  the  data  show 
substantial  scatter,  individual  points  are  accurately  determined  (see  above).  In  all  cases, 
the  observed  speeds  were  in  excess  of  the  shear  wave  speed  of  Homalite,  i.e.,  greater  than 
60%  of  the  dilatational  wave  speed.  In  general,  the  trend  is  toward  decreasing  sliding 
speed  with  increasing  static  pre-load.  However  there  is  no  clear  functional  dependence. 
Insight  provided  by  modeling  (see  below)  suggests  that  the  bands  identified  as  dominant 
in  collating  the  data  of  Fig.  6  may  have  corresponded  in  some  cases  to  the  onset  of 
interfacial  sliding,  but  in  other  cases  to  the  beginning  of  a  domain  of  increasing 
interfacial  shear  stress. 

Data  for  all  of  the  experiments  conducted  at  0.5  MPa  and  varying  impact  speed  lie  along 
the  ordinate  (the  axis  of  sliding  speed)  in  Fig.  6.  There  was  no  systematic  variation  of 
sliding  speed  with  impact  speed,  with  all  the  measured  sliding  speeds  being 
approximately  the  dilatational  wave  speed  of  Homalite. 


Sliding  Stresses  and  Loading  Rate:  The  photoelastic  fringes  can  be  equated  to  the  local 
difference  in  principal  stresses: 

cj1-cj2=F(JN/W  (1) 

where  Fa  is  the  stress-optic  coefficient  (22.6  kN/m  for  Homalite),  N  is  the  fringe  order, 
and  W  is  the  specimen  width  (in  the  through-thickness  direction,  X3).  The  stress  along  the 
center-line  of  the  Homalite  can  be  determined  from  Eq.  (1)  alone  if  the  following 
conditions  are  met:  1)  the  principal  stresses  along  the  center- line  act  in  the  x\  and  X2 
directions  and;  2)  plane  stress  conditions  exist;  and  3)  strains  S22  in  the  X2  direction  are 
negligible.  Then  Eq.  (1)  reduces  to 

a-i  (1  -v)  =  FaN/h  (2) 

where  v  is  Poisson’s  ratio  for  the  Homalite;  and  hence  the  stress  distribution  along  the 
center  of  the  fiber  can  be  determined.  The  condition  of  negligible  strain  component,  £22, 
will  not  be  met  in  the  presence  of  large  compressive  pre-stress.  Nevertheless,  Eq.  (2) 
proves  useful,  for  the  following  reasons.  First,  numerical  simulations  (see  below)  show 
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that,  in  the  absence  of  compressive  pre-stress,  dynamic  contributions  to  £n  are  very 
small;  and  the  stress,  at,  estimated  from  experiments  using  Eq.  (2)  and  that  calculated  by 
the  simulations  are  very  close.  Second,  when  the  compressive  pre-stress  is  large,  the 
additional  contribution  to  s-n  will  be  time-invariant.  Therefore,  Eq.  (2)  can  be  used  to 
estimate  the  applied  loading  rate,  dcr\ldx\,  in  steady-state  conditions  from  experimental 
data,  since  this  involves  taking  a  derivative. 

The  variation  of  stress  with  position  is  found  to  be  approximately  linear,  yielding  a 
constant  value  of  do\ldx\  for  each  test.  (Thus  experiments  and  simulations  (see  below) 
confirm  the  rapid  attainment  of  steady-state  conditions  in  the  axial  stress  and  the  proper 
application  of  Eq.  (2).)  This  steady-state  value  can  be  divided  by  the  dilatational  wave 
speed  in  Homalite  to  give  a  loading  rate  per  unit  time.  The  loading  rate  rises  in 
proportion  to  the  impact  speed  and  decreases  linearly  with  static  pre-load  (Fig.  7).  When 
the  coupling  of  the  fiber  and  the  matrix  is  increased,  the  input  momentum  is  transferred 
more  rapidly  to  the  matrix,  resulting  in  slower  loading  of  the  fiber. 

3.  Numerical  Modeling 

A  computational  simulation  of  the  test  was  set  up  in  a  commercial  finite  element  package 
(ABAQUS  1 ),  with  cohesive  elements  to  represent  interface  friction.  Symmetrical 
boundary  conditions  ( ut,  =  0;  tn  =  0)  were  imposed  along  the  center  of  the  Homalite  piece 
(x2  =  0)  and  the  top  surface  of  the  steel  piece.  A  state  of  plane  stress  was  assumed  in  the 
(jci,  X2)  plane.  Semi-infinite  elements  were  used  to  prevent  elastic  waves  bouncing  back 
from  the  right  end  of  the  model  (Fig.  8).  The  cohesive  elements  introduce  a  relationship 
between  the  shear  tractions,  Tint,  and  the  sliding  displacement  rate,  [ul  ] ,  at  the  interface. 

The  law  used  (Fig.  9)  is  intended  to  represent  friction  that  is  uniform  and  constant  in 
magnitude,  but  an  initial  linear  part  is  included  to  avoid  numerical  difficulties  where  the 
sense  of  the  motion  changes.  A  dynamic  loading  history  (Fig.  10),  corresponding  to  that 
caused  by  the  experimental  impactor,  was  imposed  as  a  prescribed  stress,  acting 

uniformly  over  the  interval  -h  <  X2  <  h  at  the  left  end  of  the  Homalite.  The  simulation 
described  in  detail  in  the  following  had  a  loading  rate,  dcrnldt  =  31.5  MPa/ps,  a  rise  time 


1  ABAQUS,  Inc.,  Pawtucket,  RI  02860 
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of  6.0  |lis,  and  linear  unloading  over  a  further  12.0  (as.  The  loading  rate  and  maximum 
stress  correspond  to  those  deduced  from  measurements  of  stress  distributions  at  the  end 
of  the  Homalite  for  the  test  with  an  impactor  velocity  of  39  m/s  and  a  compressive  pre¬ 
stress,  cr 22  =50  MPa.  The  unloading  rate  was  set  to  be  half  that  of  the  loading  rate  in 

magnitude,  which  also  corresponded  approximately  to  the  measured  rate.  Thus  the 
boundary  stress  was  not  predicted,  but  fitted  to  experiments;  but  all  other  characteristics 
of  the  deformation,  including  front  velocities  and  the  general  distribution  of  stresses,  etc., 
were  predicted  given  only  this  fitted  condition.  The  interfacial  frictional  stress  was 
assigned  the  magnitude  Zf  =  40  MPa.  This  would  correspond  to  a  friction  coefficient  of 
0.8  for  the  same  test,  if  friction  had  obeyed  Coulomb’s  law  in  the  test;  but  this,  of  course, 
may  not  have  been  the  case.  The  magnitude  of  Zf  was  held  constant  in  the  simulations. 
Since  the  loading  rate  depends  only  weakly  on  the  compressive  pre-stress  (Fig.  7),  the 
simulation  could  also  be  interpreted  as  representing  cases  of  higher  compressive  pre¬ 
stress,  if  a  lower  value  of  the  coefficient  was  believed  to  be  more  appropriate.  (The 
calculation  would  be  unchanged  because  the  coefficient  of  friction  does  not  enter 
explicitly  into  the  model.) 

3.1  General  Character  of  the  Predicted  Wave  Deformation 

Figure  1 1  shows  predicted  contours  of  the  principal  stress  difference  <7|-cr2  and  the  shear 
stress  r,2  in  the  Homalite  at  time  t  =  9.822  ps.  These  plots  and  plots  of  stress  variations 
along  single  lines  (see  below)  distinguish  four  zones  in  the  wave  propagation.  Right¬ 
most  is  the  so-called  head  wave  zone,  where  the  interface  friction  stress  is  opposite  in 
sense  to  that  in  the  trailing  zones.  (This  change  of  sign  is  not  evident  in  Fig.  1  lb,  because 
of  the  small  magnitudes  of  the  stresses  involved.)  In  the  head  wave  zone,  the  Homalite  is 
loaded  by  the  steel  and  the  wave  speed  exceeds  that  of  Homalite,  Cf,  by  a  factor  of  more 
than  2  (the  factor  being  difficult  to  pinpoint  numerically)  but  is  less  than  that  of  the  steel, 
Cm. 

Following  the  head  wave  zone  is  a  linearly  increasing  shear  stress  zone,  where  the 
interface  (friction)  stress  increases  approximately  linearly  with  position  up  to  the  limiting 
value,  Zf,  given  by  the  cohesive  law.  This  zone  derives  from  the  finite-sloped  transition 
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from  negative  to  positive  shear  stress  used  in  the  cohesive  model  to  avoid  numerical 
problems  (Fig.  9).  However,  domains  with  this  character  were  also  observed  in  the 
model  experiments,  suggesting  that  the  law  of  Fig.  9  may  fortuitously  represent,  at  least 
qualitatively,  the  true  constitutive  behaviour  of  the  interface  under  rising  shear  stress. 
The  next  zone  is  the  constant  friction  zone,  which  is  characterized  by  a  suddenly  and 
greatly  reduced  density  of  the  contours  in  o\-a2}  while  the  shear  stress  is  constant  at  the 
value,  Zf- 

The  last  zone  is  associated  with  the  unloading  part  of  the  loading  history.  The  contour 
lines  near  the  interface  become  chaotic,  which  is  also  seen  in  the  experiments  (Figs.  3  and 
4).  In  the  simulation,  alternating  contact  and  separation  zones  exist  along  the  interface  in 
this  region.  The  experimental  contours  show  similar  variations,  suggesting  the  same 
contact  behaviour  (although  in  the  experiments  no  direct  confirmation  of  interface 
displacements  is  possible). 


4.  Stress  Distributions  and  Front  Velocities 

Comparison  of  the  measured  and  predicted  stress  contours  has  shown  that,  even  though 
the  numerical  model  is  based  on  an  ad  hoc  and  probably  incorrect  friction  law,  most 
aspects  of  the  experiments  appear  to  be  reproduced.  Further  insight  into  the  physics  of 
the  deformation  is  available  from  other  calculated  and  measured  characteristics. 

Figures  12a  -  12c  show  stress  profiles  at  three  different  instants  during  the  simulations  of 
Fig.  11.  Figure  12a  shows  the  axial  stress  along  the  center  line  of  the  Homalite,  crn(xi, 
0),  along  with  a  measurement  from  the  experiment  executed  with  impactor  velocity,  vo  = 

39  m/s  and  compressive  pre-stress,  cr^  =  50  MPa.  The  agreement  between  the 

simulation  and  the  experiment  is  quite  close.  (See  also  the  discussion  following  Eq.  (2).) 
The  furthest  disturbance  evident  in  this  plot  corresponds  to  the  transition  in  Fig.  11  from 
the  head  wave  to  the  domain  of  linearly  increasing  friction.  The  peak  axial  stress  along 
the  center  line  decreases  as  the  wave  propagates,  which  is  attributed  to  the  dissipative 
frictional  sliding  process  along  the  interface.  Invariance  of  the  stress  profile  confirms  the 
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attainment  of  steady  state  conditions  by  time  6.2  psec,  shortly  after  peak  load  is  attained 
at  6.0  psec. 

Figure  12b  shows  the  interface  shear  stress  at  the  same  three  instants.  The  shear  stress 
increases  approximately  linearly  until  the  constant  friction  stress  assigned  in  the  cohesive 
law  is  reached.  It  remains  almost  constant  through  the  “constant  friction”  zone  of  Fig. 
1 1 ,  although  small  oscillations  are  evident  around  the  prescribed  value,  Tf.  The  unloading 
zone  is  dominated  by  high  frequency,  severe  fluctuations  in  r;n t. 

Figure  12c  shows  the  axial  stress,  <J\\(x\,  h ),  along  the  interface.  A  knee  behind  the 
furthest  visible  disturbance  corresponds  to  the  transition  from  linearly  increasing  to 
constant  interfacial  stress.  In  the  region  ahead  of  the  knee,  the  loading  rate  of  the 
material  is  remarkably  small. 

Figure  12d  illustrates  further  the  nature  of  the  interface  fields  in  the  chaotic  unloading 
zone.  Large  spikes  in  the  fiber  sliding  displacement  periodically  reach  down  to  the 
sliding  displacement  of  the  matrix,  indicating  locations  of  fiber/matrix  stick.  Reverse  slip 
is  also  possible  for  brief  intervals  at  these  locations,  but  is  difficult  to  resolve  in  a 
numerical  simulation,  due  to  computational  noise. 

The  plots  of  Fig.  12  provide  reasonably  accurate  information  for  evaluating  the  velocities 
of  the  various  fronts.  Fig.  13  shows  histories  of  the  locations  of  fronts  as  functions  of 
time,  obtained  by  interpolating  measurements  of  the  locations  of  the  corresponding 
features  in  Fig.  12.  The  leading  edge  of  the  head  wave  is  not  shown,  since  it  could  not  be 
accurately  discerned. 

The  velocity,  V2,  of  the  transition  to  the  linearly  increasing  friction  zone  is  evaluated  from 
both  the  earliest  significant  axial  stress  along  the  center  line  and  the  earliest  significant 
shear  stress  along  the  interface.  The  former  gives,  V2  =  1 ,04Cf;  while  the  latter  yields  a 
slightly  lower  value  V2  =  0.98Cf.  To  within  reasonable  uncertainty,  neither  is 
significantly  different  from  the  wave  speed  in  Homalite. 

The  front  of  the  constant  friction  zone  can  only  be  determined  accurately  from  the 
interface  stresses,  either  the  knee  in  the  axial  stress  or  the  onset  of  the  shear  plateau  (Fig. 
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12b).  The  inferred  wave  velocity  is  V3  =  0.87  ±  0.02  Cf,  the  error  reflecting  the 
difference  between  the  two  sources  and  the  axial  stress  implying  the  higher  value. 

The  front  associated  with  the  onset  of  unloading  is  best  identified  by  the  peak  in  the  axial 
stress  along  the  center  line  of  the  Homalite  (Fig.  12c).  The  velocity  takes  the  value  V4  = 
0.98  ±  .02  Cf,  once  again  close  to  the  wave  speed  in  the  Homalite. 

The  velocity  estimates  can  be  approximately  confirmed  by  analysing  the  bands  or  kinks 
in  the  predicted  stress  contours,  which  propagate  from  locations  on  the  interface  at  which 
a  transition  (non-smooth  behaviour)  is  present  in  the  displacement  fields.  Two  bands  are 
traced  by  dashed  lines  in  Fig.  11,  emanating  approximately  from  the  fronts  associated 
with  the  onset  of  linearly  increasing  friction  and  constant  friction  (sliding).  The  angles 
subtended  by  the  lines  to  the  interface,  <25,  i  =  2  and  3,  are  given  by  sin  a,  =  Cs/Fi,  where 
Cs  is  the  shear  wave  speed  in  the  Homalite.  While  difficulty  arises  in  locating  the  bands 
on  some  contours,  the  fitted  values  sinai  «  0.45  and  sinas  «  0.6  are  probably  correct  to 
within  20%.  These  values  yield  F?  =  (2.2  ±  0.4)Cs  and  Cs  =  (1.6  ±  0.3)Cs,  in  reasonable 
agreement  with  the  more  accurate  values  deduced  above  from  Fig.  12.  The  boundary  of 
the  domain  of  chaotic  fields  is  also  marked  by  a  dashed  line  in  Fig.  11,  which  emanates 
from  the  unloading  front.  The  angle  subtended  by  this  line,  a4  «  arcsin  0.35,  yields  FChaos 
=  (3.0  ±  0.3)Cs.  This  high  implied  velocity  is  consistent  with  the  observation  in  the 
simulations  that  the  onset  of  chaos  is  delayed  somewhat  after  the  beginning  of  unloading, 
but  the  chaotic  zone  tends  subsequently  to  catch  up  with  the  unloading  front,  and 
therefore  propagates  at  a  velocity  exceeding  that  of  the  unloading  front,  or  approximately 
Cf. 

The  different  velocities  predicted  for  the  onset  of  sliding  and  the  beginning  of  the  domain 
of  increasing  interface  shear  help  understand  the  distribution  of  the  data  in  Fig.  6.  Some 
uncertainty  arose  in  identifying  bands  in  fringe  images  (as  described  above),  so  that  the 
nature  of  the  interfacial  phenomenon  corresponding  to  a  particular  band  was  not  always 
clear.  Given  the  velocities  predicted  by  the  model,  those  data  in  Fig.  6  falling  close  to  Cf 
would  be  associated  with  the  beginning  of  the  domain  of  rising  interfacial  shear  stress, 
while  other  data  would  be  associated  with  the  sliding  front. 


17 


Since  the  front  of  the  constant  friction  zone  propagates  at  slightly  different  velocities  at 
the  interface  and  along  the  center  of  the  Homalite,  the  configuration  of  contours  must 
change  with  time.  Other  contour  plots  at  different  times  (not  shown)  show  that  the  nearly 
parallel,  vertical  contours  in  the  “linearly  increasing  friction”  and  “constant  friction” 
zones  of  Fig.  1 1  do  move  apart  more  rapidly  at  the  center  of  the  Homalite  than  at  the 
interfaces,  which  causes  them  to  bow  out  towards  the  head  wave  with  increasing  time. 
The  same  fringe  divergence  was  seen  (via  dynamic  photoelasticity)  in  the  model 
experiments.  In  both  the  simulations  and  the  experiments,  the  divergence  is  limited  in 
magnitude  and  the  stress  contours  remain  at  angles  less  than  approximately  10°  from 
vertical  during  passage  of  the  deformation  along  the  length  of  the  test  piece. 

Other  changes  with  time  are:  1)  The  zone  of  chaotic  contours  grows  into  the  interior  of 
the  Homalite  with  time.  In  fact,  other  contour  plots  from  the  same  simulations  as 
reported  in  Fig.  1 1  reveal  that  chaotic  behaviour  does  not  appear  at  all  until  the  elapse  of 
approximately  8  ps.  2)  The  constant  friction  stress  zone  shrinks  with  time,  since  V4  >  V3, 
indicating  that  the  unloading  wave  propagates  faster  than  waves  ahead  of  it. 


5.0  Results  from  Shear  Lag  Models  based  on  Lame  Fields 

For  the  static  problem  of  Fig.  1  and  related  thermal  and  mechanical  problems,  many 
useful  results  have  been  obtained  using  shear  lag  analysis  and  the  assumption  that  the 
fiber  deformation  fields  are  separable  in  x  and  the  transverse  or  radial  variables;  and  that 
fields  are  simple,  known  functions  of  the  latter  (Cox,  1990;  Cox  et  al.,  1990;  Hutchinson 
and  Jensen,  1990;  Marshall  et  al.,  1985;  McCartney,  1987).  Analysis  shows  that  these 
approximations  give  quite  accurate  predictions  of  the  end  load  vs.  the  end  displacement 
as  long  as  the  slip  distance  is  long  compared  to  the  fiber  diameter  (or  fiber  width  in  a 
plane  problem)  (Hutchinson  and  Jensen,  1990).  The  fiber  fields  are  well  represented;  and 
also  the  matrix  fields,  provided  the  volume  fraction  of  fibers  is  not  too  small.  Since  these 
results  in  the  static  case  are  simple  and  have  enlightened  many  aspects  of  the  mechanics 
and  engineering  principles  of  designing  fiber-reinforced  materials,  the  question  of  their 
utility  in  the  dynamic  case  is  of  interest. 

For  dynamic  end  loading  that  increases  linearly  in  time  and  friction  that  is  uniform  and 
constant,  shear  lag  analysis  predicts  that  the  deformation  will  propagate  in  two  or  three 
domains,  depending  on  the  value  taken  by  the  following  three  dimensionless  parameters 
(Sridhar  et  al.,  2003): 

k=T'C‘  C  =  A_  &-  (4) 

h  !  dt  c„  r  (1  -/K, 
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where  /  is  the  fiber  volume  fraction.  The  specimen  of  Fig.  2  and  the  simulations  have 
lateral  boundary  conditions  that  are  not  periodic  and  therefore  are  not  in  strict 
correspondence  with  the  shear  lag  model  of  Sridhar  et  al.,  but  correspondence  can  still  be 
sought,  using  the  dimensions  of  the  specimen,  by  setting /=  8.25/50.8  =  0.162.  The 
domains  predicted  by  shear  lag  modeling  consist  of  different  interfacial  conditions:  slip, 
stick,  and  reverse  slip.  For  the  case  of  the  simulations  reported  in  the  last  section,  k  = 
0.317,  C  =  0.378,  and  (f>  =  0.0043.  Shear  lag  modeling  predicts,  for  these  values,  that  the 
deformation  during  the  linearly  rising  part  of  the  loading  will  comprise  a  stick  zone  and  a 
slip  zone  (see  Appendix  A).  Furthermore,  the  interface  friction  stress  in  the  stick  zone  is 
negative  and  has  a  very  small  magnitude,  0.036  MPa  (see  Eq.  (A.6)).  The  sign  is 
negative  because  the  matrix  (steel),  having  the  higher  wave  speed,  is  pulling  the  fiber 
along  in  this  zone.  The  magnitude  is  small  because  the  volume  fraction  of  the  fiber  is 
small,  so  that  the  stresses  in  the  matrix  remain  small,  the  spatial  derivative  remains  small, 
and  therefore  the  friction  stress  remains  small.  For  the  case  studied,  the  shear  lag  model 
predicts  (Appendix  A)  that  the  head  wave  front  will  propagate  at  V\  =  2.6  Cf. 

A  further  prediction  of  the  shear  lag  model  is  that,  if  the  fiber  volume  fraction  rises,  then 
the  interface  shear  stress  in  the  head  wave  zone  will  exceed  the  friction  stress  in 
magnitude  and  reverse  slip  will  occur  (consider  Fig.  A.l  as  the  parameter  (f>  rises).  The 
head  wave  zone  would  no  longer  be  a  stick  zone. 

Figure  1 1  broadly  confirms  the  predictions  concerning  the  head  wave  zone.  The  zone  is  a 
stick  zone  (no  slip),  the  interface  shear  stress  is  negative  and  very  small  in  magnitude, 
and  the  front  velocity,  while  not  well  determined,  is  consistent  with  the  value  2.6  Cf. 

The  zone  of  linearly  increasing  friction  in  Fig.  11  has  no  analogue  in  the  shear  lag  results 
of  Sridhar  et  al.,  because  in  that  particular  modeling,  the  friction  stress  was  assumed  to 
change  sign  with  reversing  slip  direction  as  a  step  function.  The  zone  of  linearly 
increasing  friction  in  Fig.  11  arises  from  the  linear  friction  stress  regime  at  small 
displacements  in  the  law  of  Fig.  9.  The  constant  friction  zone  of  Fig.  1 1  is  equivalent  to 
the  slip  zone  in  the  shear  lag  model.  The  front  of  the  constant  friction  zone  is  predicted 
in  the  shear  lag  model  to  propagate  at  0.73  Cf.  This  is  less  than  the  velocity,  V3  =  0.87  Cf, 
computed  in  the  simulations  of  Fig.  11.  An  interesting  question  is  whether  the  two 
models  would  agree  in  the  predicted  front  velocity  for  higher  fiber  volume  fractions, 
equivalent  lateral  boundary  conditions,  and  as  the  linearly  varying  domain  in  the 
constitutive  law  of  Fig.  9  shrinks  in  extent. 

If  it  is  assumed  that  the  friction  stress  value  used  in  the  shear  lag  model  would  vary  with 
compressive  pre-stress  according  to  Coulomb’s  law,  then  a  prediction  of  the  variation  of 
front  velocities  with  compressive  pre-stress  can  be  made.  Two  such  predictions  for  the 
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velocity  of  the  sliding  front,  V3,  have  been  superimposed  on  Fig.  6.  The  rates  of  change 
of  the  friction  stress  correspond  to  assuming  a  coefficient  of  friction  of  0. 1  or  0.4.  The 
trends  of  the  curves  are  not  inconsistent  with  the  weakly  defined  declining  trend  in  the 
data,  especially  when  those  data  lying  at  Cf,  which  are  believed  to  correspond  to  the  front 
of  increasing  interfacial  shear  stress,  are  excluded. 

6.0  Possible  Chaotic  Behaviour 

The  noisy  behaviour  appearing  in  the  unloading  zone  in  the  simulations  does  not  appear 
to  a  numerical  artefact,  but  rather  a  consequence  of  the  modeling  assumptions  and  the 
physics  of  the  problem.  If  it  were  due  to  Gibbs’  phenomenon,  i.e.,  numerical  ringing  at  a 
near-discontinuity  in  the  solution,  it  would  not  be  confined  to  the  vicinity  of  the  interface 
(note,  for  example,  the  smoothness  of  the  axial  stress  at  the  center-line  of  the  Homalite, 
Fig.  12a);  nor  to  the  unloading  zone,  since  non-smoothness  of  equal  strength  exists  in  the 
fields  at  other  fronts  (zone  boundaries). 

The  noisy  behaviour  is  not  associated  with  Coulomb’s  law,  since  this  was  not  the 
constitutive  assumption  of  the  model.  In  fact,  the  friction  law  assumed  (Fig.  9),  which 
contains  no  relation  between  the  friction  stress  and  the  normal  pressure,  is  simpler  than 
that  used  in  most  studies  of  either  stability  or  ill-posedness  in  problems  of  slip  between 
half-spaces.  For  the  same  reason,  the  predicted  phenomenon  cannot  be  directly  related  to 
Poisson’s  effect,  even  though  Poisson’s  ratio  was  nonzero  in  the  simulations. 

The  fact  that  the  noisy  behaviour  is  confined  to  the  unloading  zone,  while  the 
deformation  propagates  smoothly  at  a  range  of  speeds  in  other  zones,  has  no  obvious 
analogue  in  any  prior  studies  of  instability.  The  closest  to  a  hint  of  similar  behaviour  is  in 
incomplete  and  unpublished  results  for  the  shear  lag  modeling  described  by  (Sridhar  et 
al.,  2003).  Attempts  to  trace  solutions  for  linear  unloading  that  follows  linear  loading 
lead  to  apparently  complex  patterns  of  the  birth  and  death  of  slip,  stick,  and  reverse-slip 
zones,  as  the  unloading  wave  overtakes  deformation  of  earlier  origin.  The  possibility  is 
suggested  (but  not  yet  thoroughly  researched)  that  arbitrarily  small  changes  in  the  loading 
history  causing  finite  changes  in  the  subsequent  pattern  of  slip,  stick,  and  reverse-slip 
zones.  In  the  shear  lag  problem,  the  only  possible  source  of  instability  or  complexity  is 
the  strong  nonlinearity  associated  with  the  change  in  sign  of  the  friction  stress  when  the 
direction  of  interfacial  slip  changes.  While  the  analysis  is  incomplete,  the  suggestion  of 
chaotic  behaviour  in  the  shear  lag  model  provokes  the  speculation  that  the  noisy 
behaviour  seen  during  unloading  in  the  experiments  and  the  numerical  simulations  is  also 
an  example  of  chaos. 

Shear  lag  modeling  also  highlights  the  role,  in  the  possible  generation  of  chaotic 
behaviour,  of  the  coupling  of  wave  motions  between  the  fiber  and  the  matrix.  For  a  fiber 
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in  a  rigid  matrix,  where  the  fiber  deformation  is  modified  by  interfacial  friction  tractions 
that  depend  only  on  its  own  motion,  stable  and  well-behaved  solutions  can  be  mapped  out 
systematically  for  quite  complex  loading  and  unloading  histories,  e.g.,  by  the  method  of 
characteristics  (Nikitin  and  Tyurekhodgaev,  1990).  This  method  is  not  useful  when  the 
fiber  motion  is  coupled  to  that  of  a  matrix. 

Numerical  solution  of  the  shear  lag  problem,  e.g.,  numerical  solving  of  difference 
equations  subject  to  the  constraint  that  displacement  fields  are  Lame-like,  is  not 
enlightening,  because  the  rapid  births  and  deaths  of  different  zones  are  quickly  obscured 
in  numerical  noise  at  the  fronts. 

The  similarity  of  the  predicted  domains  of  chaotic  behaviour  and  those  measured  in  the 
model  experiments  is  very  eye-catching  (Figs.  4  and  11),  but  the  conditions  assumed  for 
the  interface  are  unlikely  to  be  the  same  in  the  two  cases.  In  particular,  the  friction  law  of 
Fig.  9  is  an  idealization  that  is  likely  to  be  an  oversimplification  of  the  true  friction 
behaviour  in  the  experiments  (although  the  experimental  law  is  not  known  in  detail). 
This  suggests  that  the  putative  existence  of  chaotic  behaviour  is  not  especially  sensitive 
to  details  of  the  friction  law,  but  could  be  a  consequence  of  unloading  with  any  friction 
law  that  possesses  the  strong  nonlinearity  of  sign  reversal. 


7.0  Conclusions 

Comparing  the  characteristics  of  model  planar  experiments,  numerical  simulations,  and 
shear  lag  modeling  leads  to  the  following  inferences  about  the  push-in  problem  of  interest 
(Fig.  1),  for  a  bi-linear  sequence  of  loading  and  unloading. 

1.  The  broad  characteristics  of  the  deformation  in  the  fiber,  exclusive  of  shock 
waves,  can  be  predicted  at  least  qualitatively  by  simple  shear  lag  models,  in  which  the 
displacements  are  reduced  to  Lame  fields.  The  shear  lag  solutions  have  the  advantage  of 
proposing  clearly  defined  zones  of  different  slip  behaviour,  which  while  present  in  both 
the  experiments  and  the  numerical  simulations,  are  not  so  easily  identified  in  them 
without  prior  knowledge.  For  engineering  studies  of  dynamic  deformation  in  composites 
reinforced  by  fibers,  rods,  or  stitches,  shear  lag  results  will  convey  at  least  the  trends  of 
dynamic  response  at  the  macroscopic  level.  For  example,  the  dynamic  fiber  end 
displacement  is  likely  to  be  described  quite  well  by  shear  lag  theory,  since  those  fine 
details  of  the  push-in  (or  pullout)  phenomenon  that  shear  lag  theory  misses  are  unlikely  to 
strongly  influence  an  quantity  such  as  the  end  displacement,  which  is  derived  from 
integrated  strains. 
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2.  For  the  given  geometrical  conditions  and  loading  rate,  shear  lag  modeling 
correctly  predicts  the  system  of  stick  followed  by  slip  found  in  the  numerical  simulations 
and  model  experiments,  even  though  the  details  of  the  assumed  or  actual  friction  laws  are 
likely  to  be  different  in  all  three.  One  infers  that  the  pattern  of  zones  is  mainly 
determined  by  the  interplay  of  stress  waves  of  different  speeds  in  the  fiber  and  the  matrix 
and  is  fairly  insensitive  to  the  friction  constitutive  law. 

3.  The  speeds  of  fronts  (or  the  boundaries  between  zones)  are  fairly  well  predicted 
by  shear  lag  modeling,  but  more  consistent  between  the  numerical  simulations  and  the 
model  experiments.  Shear  lag  modeling  could  be  expected  to  be  a  better  approximation 
as  the  volume  fraction  of  the  fibers  rises. 

4.  As  far  as  shear  lag  theory  is  accurate,  results  found  here  for  plane  geometry 
should  be  equally  applicable  to  axisymmetric  conditions,  which  are  a  useful 
approximation  to  composites  of  cylindrical  fibers. 

5.  The  literature  is  replete  with  analyses  of  frictional  sliding  between  half- spaces 
under  uniform  far-field  velocity  conditions.  For  that  problem,  the  presence  of  ill- 
posedness,  instability,  pulse  and  crack-like  slip  systems,  etc.,  depends  strongly  on  the 
nature  of  the  friction  law  (see  Introduction).  The  present  study  considers  non-uniform 
(time  dependent)  loading,  which  yields  results  that  are  not  accounted  for  by  the  prior 
literature.  Complex  behaviour,  which  may  be  chaotic,  arises  during  unloading,  but  not  at 
other  stages  of  the  deformation,  in  both  the  experiments  and  the  numerical  solutions. 
This  behaviour  appears  to  be  a  consequence  of  the  strong  nonlinearity  associated  with  the 
sign  reversal  in  the  frictional  tractions  with  slip  direction  and  is  not  sensitive  to  the  details 
of  the  friction  constitutive  law. 

6.  The  results  that  some  characteristics  of  engineering  relevance,  especially  the  fiber 
end  displacement  and  the  pattern  and  velocity  of  sliding  fronts,  are  insensitive  to  details 
of  the  friction  law  (being  so  similar  in  the  experiments,  numerical  modeling,  and  shear 
lag  analysis)  also  contrasts  with  the  literature  on  shear  crack  propagation.  One  would 
expect  that  as  the  frictional  traction  rises  in  magnitude,  or  a  large  debond  energy  is 
included  in  the  problem,  the  fiber  push-in  problem  would  begin  to  manifest  the  complex 
dependence  seen  for  shear  cracks.  At  the  same  time,  the  conditions  required  for  shear  lag 
analysis  to  be  accurate  would  be  lost. 

7.  Changing  the  compressive  contact  pressure  in  the  experiments  has  quantitative 
but  not  qualitative  effect  on  the  deformation  patterns. 
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Appendix  —  results  from  shear  lag  analysis 

The  following  results  are  reproduced  from  (Sridhar  et  al.,  2003)  and  apply  to  the 
fiber/matrix  geometry  of  Fig.  1  when  friction  is  uniform  and  constant  and  the  fiber  is 
subjected  to  a  linearly  increasing  end  load.  For  this  problem,  three  cases  arise  in  the 
solutions,  viz.,  cases  where  domains  of  slip  occur,  cases  where  domains  of  slip  and  stick 
occur,  and  cases  where  domains  of  slip  and  reverse  slip  occur.  Figure  A.l  shows  how 
which  of  these  cases  prevails  depends  on  the  loading  rate  and  material  parameters  of  the 
problem.  The  parameter 


k  =  Tf  °i  j  d(To 
h  /  dt 


(A.1) 


represents  the  loading  rate,  with  Cf  the  longitudinal  wave  speed  in  the  fiber  (Homalite),  Tf 
the  constant  friction  stress,  and  do-^ldt  the  rate  of  increase  of  the  stress  at  the  fiber  end. 
With  the  subscripts  f  and  m  referring  to  the  fiber  and  matrix,  respectively,  the  parameters 
C  and  (f>  are  defined  by 
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with  Ex  =  El  (l  -  //£  )/((l  +  //;  )(l  -  2 //j )) ,  i  =  m  or  f;  and  Eh  //,,  and  p,,  i  =  m  or  f,  denoting 
Young’s  modulus,  Poisson’s  ratio,  and  the  density. 

The  case  of  the  test  shown  in  Figs.  3,  10,  and  11,  i.e.,  for  38  m/s  impact  speed,  falls  into 
the  regime  of  stick/slip  in  the  map  of  Fig.  A.l.  The  front  of  the  stick  zone  advances  at  a 
speed  V\  given  by 


n=cf 


\+(f) 

c2+</> 


(A.4) 


while  the  front  at  which  slip  begins  advances  at  a  speed  V3  given  by  the  root  of  a  cubic 
equation  in  773  =  Vj/cf. 

(773  +2*73  -1)(1  +  C2773771)  +  2  k  p(l  +  737l)%  =°  (A-5) 
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Analysis  shows  that  73  has  only  one  real  positive  root  which  always  lies  in  (0,1),  whereas 
7 1  can  clearly  exceed  unity.  The  interfacial  friction  stress  in  the  stick  zone  is: 

.  (y|  +2A>?3 -1)(1-C2) 

2  *  (1  -  c2 )  +  (1  -  nl  )V(C2  +^)(1+(CT) 

where  73  has  been  obtained  by  solving  Eqn.  (A.5).  This  expression  for  the  friction  stress 
always  satisfies  |  r  |  <  14 


(A.6) 
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periodic  boundary  conditions  (plane  case) 

or  type  I  or  type  II  boundary  conditions  (axisymmetric  case) 


dynamic 
end  load 


Figure  1.  Schematic  of  a  general  problem  of  fiber/matrix  debonding,  under  dynamic 
end  loading,  with  possible  crack  wake  friction. 


Figure  2:  Planar  test  specimen  under  static 
pre-load  applied  to  the  top  and  bottom  plates 
of  the  fixture. 
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T  constant 


Figure  3.  Representative  image  from  a  dynamic  sliding  experiment  (38  m/s  impact 
speed,  50  MPa  static  pre-load).  (The  constancy  of  the  interface  shear  stress  in  the 
region  labeled  “t  constant”  is  consistent  with  the  fringe  pattern,  but  not  well  resolved 
experimentally  above  the  noise.  It  is  marked  so  because  of  modeling  results.  The 
labels  of  the  other  domains  are  clearly  implied  by  the  nature  fringe  patterns.) 


Figure  4.  Fringe  patterns  obtained  at  low  and  high  values  of  the  static  pre-load  (as 
marked). 
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time  (jus) 

Figure  5:  A  plot  of  the  location  of 
sliding  onset  vs.  time  for  a  specimen 
impacted  at  40  m/s  with  a  static  pre¬ 
load  of  12.5  MPa. 


confinement  pressure  (MPa) 

Figure  6:  A  plot  of  the  variation  of  the 
speed  of  the  sliding  front  with  static  pre¬ 
load. 
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Figure  7:  Variation  of  loading  rate  with 
magnitude  of  static  compressive  pre¬ 
load. 
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Figure  8:  Mesh  for  numerical  simulations. 


Figure  9:  Interface  friction  law  introduced 
via  cohesive  elements. 


time,  t  (fis) 


Figure  10.  Loading  history  in 
simulation. 
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X1  (mm) 

Figure  11:  Calculated  contour  plots  for  (JX-<J2  and  r12  in 
the  Homalite  at  t  =  9.822  ps.  Arrows  indicate  direction  of 
increasing  magnitude  of  stress. 
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(a)  (b) 


x  (mm) 
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(c)  (d) 

Figure  12.  Stress  and  displacement  profiles  from  numerical  simulations,  (a)  axial  stress, 
an,  along  the  center  line  of  Homalite  (heavy  curve  is  experimental  data);  (b)  shear  stress, 
X12;  (c)  axial  stress,  an,  along  the  interface;  and  (d)  axial  particle  velocity  in  the  Homalite 
(fiber)  and  steel  (matrix)  along  the  interface. 
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axial  front  along  center 


Figure  13.  Front  motions. 


(P  - ► 

Fig.  A.l.  Domains  in  which  different  configurations  of  slip,  stick,  and  reverse 
slip  are  predicted  by  shear  lag  analysis  for  a  fiber  push-in  problem  with 
constant  loading  rate 
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